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Abstract 

As demonstrated by Slepian et. al. in a sequence of classical papers (see [33], [33], [17] . 
[35] . [36]). prolate spheroidal wave functions (PSWFs) provide a natural and efficient tool for 
computing with bandlimited functions defined on an interval. As a result, PSWFs are becoming 
increasing popular in various areas in which such function occur - this includes physics (e.g. 
wave phenomena, fluid dynamics), engineering (e.g. signal processing, filter design), etc. 

To use PSWFs as a computational tool, one needs fast and accurate numerical algorithms for 
the evaluation of PSWFs and related quantities, as well as for the construction of quadratures, 
interpolation formulas, etc. For the last half a century, substantial progress has been made in 
design of such algorithms - this includes both classical results (see e.g. [3]) as well as more 
recent developments (see e.g. [38]1. 

The complexity of many of the existing algorithms, however, is at least quadratic in the 
band limit c. For example, the evaluation of the nth eigenvalue of the prolate integral operator 
requires at least 0(c 2 ) operations (see e.g. [3S]); also, the construction of accurate quadrature 
rules for the integration of bandlimited functions of band limit c requires 0(c 3 ) operations (see 
e.g. [6]). Therefore, while the existing algorithms are quite satisfactory for moderate values of 
c (e.g. c < 10 3 ), they tend to be relatively slow when c is large (e.g. c > f0 4 ). 

In this paper, we describe several numerical algorithms for the evaluation of PSWFs and 
related quantities, and design a class of PSWF-based quadratures for the integration of bandlim- 
ited functions. Also, we perform detailed analysis of the related properties of PSWFs. While 
the analysis is somewhat involved, the resulting numerical algorithms are quite simple and ef- 
ficient in practice. For example, the evaluation of the nth eigenvalue of the prolate integral 
operator requires 0(n + c) operations; also, the construction of accurate quadrature rules for 
the integration of bandlimited functions of band limit c requires O(c) operations. 

Our results are illustrated via several numerical experiments. 
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1 Outline 



1.1 Quadratures for Bandlimited Functions 

The principal goal of this paper is a quadrature designed for the integration of bandlimited functions 
of a specified band limit c > 0. 

A function / : R — > K is bandlimited of band limit c > 0, if there exists a function a e L 2 [—1,1] 
such that 

f(x) = j a(t) ■ e icxt dt. (1) 

In other words, the Fourier transform of a bandlimited function is compactly supported. While (UJ 
defines / for all real x, one is often interested in bandlimited functions, whose argument is confined 
to an interval, e.g. — 1 < x < 1. Such functions are encountered in physics (wave phenomena, fluid 
dynamics), engineering (signal processing), etc. (see e.g. [33], [TO], [2"0]). 
By quadrature we mean a set of nodes 

-1 < t[ n) < ■ ■< t^ < 1 (2) 

and weights 

wi n) ,...,w<r>. (3) 

If / : (— 1, 1) —> K is a bandlimited function, we use the quadrature to approximate the integral of 
/ over the interval (— 1, 1) by a finite sum; more specifically, 

/l n 
/(t)*«E w i n) /(*S n) )- (4) 

- 1 4=1 

About half a century ago it was observed that the eigenfunctions of the integral operator F c : 
L 2 [—1,1] — > L 2 [—1,1], defined via the formula 

F c [<p] (x) = [ v(t)e lcxt dt, (5) 



provide a natural tool for dealing with bandlimited functions, defined on the interval [—1,1]. More- 
over, it was observed (see [33], [T7], [35]) that the eigenfunctions of F c are precisely the prolate 
spheroidal wave functions (PSWFs) of band limit c, well known from the mathematical physics (see, 
for example, [53], [TO]). Therefore, when designing a quadrature for the integration of bandlimited 
functions of band limit c > 0, it is natural to require that this quadrature integrate several first 
PSWFs of band limit c with high accuracy. 

We formulate the principal objective of this paper in a more precise manner, as follows. 

Principal goal of this paper. Suppose that c > is a real number. For every integer n > 0, 
we define a quadrature of order n (for the integration of bandlimited functions of band limit c over 
(—1,1)) by specifying n nodes and n weights (see ([2]), ©J. Suppose also that e > 0. We require 
that, for sufficiently large n, the quadrature of order n integrate the first n PSWFs of band limit c 
up to the error e. More specifically, we find the integer M = M(c, e) such that, for every integer 
n > M and all integer m = 0, 1, . . . , n — 1, 



1 n 

1 .7 = 1 



< e, (6) 



where ip m : (—1, 1) — > K is the mth PSWF of band limit c (see Section \2.1)) . 
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Quadratures for the integration of bandlimited functions which satisfy ([BJ have already been 
discussed in the literature, for example: 

Quadrature 1. Suppose that n > is an integer. The existence and uniqueness of n nodes and 
weights, such that ([BJ holds for e — and all m — 0, 1, . . . , 2n — 1, was first observed more than 
100 years ago (see, for example, [15] . [IB] . [21] , |22j ) for all Chebyshev systems, of which PSWFs are 
a special case (see [38]). Although numerical algorithms for the design of this optimal quadrature 
were recently constructed (see [BJ, [10], [SH])) they tend to be rather expensive (require order n 3 
operations with a large proportionality constant). 

Quadrature 2. Another quadrature was suggested in [38]. The PSWF ip n has n roots t±, . . . ,t n 
in the interval (—1, 1) (see Theorem[T]in Section l!0]) ; the idea is to use these roots as the quadrature 
nodes, solve the linear system of n equations 



for the unknowns W\ , . . . , W n , and use the resulting weights and nodes to define a quadrature for 
the integration of functions of band limit 2c. This approach is justified by the generalization of 
the Euclid's division algorithm for PSWFs (see [38]), and is less expensive computationally than 
the previous one (its cost is dominated by the cost of solving the linear system ([7])). The same 
quadrature can be used to integrate functions of band limit c, since ([7]) implies that ([BJ holds with 
e = 0, for all m = 0, . . . , n — 1. 

In this paper, we describe another quadrature whose nodes are the n roots of tp n in (—1,1). 
However, its weights differ from the solution of ([7j, and can be evaluated in 0(n) operations (see 
Section rOl and Section [5] below). 

Thus, the quadratures of this paper are much faster to evaluate than those described above. 
Moreover, ([BJ ensures that their accuracy is similar to that of Quadrature 2. Also, their nodes and 
weights can be used as starting points for the scheme that computes the optimal Quadrature 1. 

In order to define the weights, to make sure that ([BJ holds and to be able to compute them 
efficiently, we need to analyze the PSWFs in a somewhat detailed manner. This analysis will be 
preceded by a heuristic explanation, which provides some intuition as well as prevents one from the 
danger of not seeing the forest for the trees (see Section [TT21 below) . Section [TT3I contains a short 
overview of the analysis. Section [5] contains mathematical and numerical preliminaries, to be used 
in the rest of the paper. In Section [3J we summarize the principal analytical results of the paper. 
Section [?] contains the corresponding theorems and proofs. Section [5] contains the description and 
analysis of the numerical algorithms for the evaluation of the quadrature and some related quantities. 
In Section [BJ we report the results of several numerical experiments. 

1.2 Intuition Behind Quadrature Weights 

We recall the following classical interpolation problem. Suppose that t±,...,t n are n distinct points 
in the interval (—1, 1). We need to find the real numbers W\, . . . , W n such that 



for all the polynomials p of degree at most n—1. In other words, the quadrature with nodes t\, . . . , t n 
and weights W\ , . . . , W n integrates all the polynomials of degree up to n — 1 exactly (see @ , (J3J , 





(8) 



©)• 
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To solve the problem, one constructs n polynomials l\, . . . , l n of degree n — 1 with the property 

lift) = I? (9) 




for every integer i, j = 1, . . . , n (see, for example. |14j). It is easy to verify that, for every j — 1, . . . , n, 
the polynomial lj is defined via the formula 

h<t)= l(f w . n{ f! (io) 
<(*i) ■ {t-tj) 

where w n is the polynomial of degree n whose roots are precisely t%, . . . , t n . The weights Wi, . . . , W n 
are defined via the formula 




for every integer j = 1, . . . , n. We observe that a single function w n is used to define all the n weights; 
also, uu n is a polynomial of degree n, and hence does not belong to the space of the polynomials of 
degree up to n — 1. 

In our case, the basis functions are the PSWFs and not the polynomials. Suppose that the roots 
ii, . . . , t n of ip n in the interval (—1, 1) are chosen to be the nodes of the quadrature. If we choose the 
weights Wi, . . . , W n such that the resulting quadrature integrates the first n PSWFs exactly, this 
will lead to the linear system (J7J, and hence to Quadrature 2 from Section [01 Instead, we define 
the weights via using ip n in the same way we used w n in (jlll) . More specifically, similar to (|10j) . for 
every integer j = 1, . . . , n, we define the function ipj : (— 1, 1) — > R via the formula 

ft* = 77771 77 — 7T- ( 12 ) 
We observe that, for every integer i, j = 1, . . . , n, 

In t -L a 

(13) 

analogous to © . Viewed as a function on the whole real line, each ipj is bandlimited with the same 
band limit c (see, for example, Theorem 1591 in Section l4.4.1l or Theorem 19.3 in |31jV On the other 
hand, <pj does not belong to the span of . . . , tp n -i (see Theorem 1551 in Section |4.4.1[) . We 

define the weights W\ , . . . , W n via the formula 

Wj = J_<piV) dt > ( 14 ) 

for j — 1,2, ... ,n. The weights Wi, . . . , W n , defined via (TT41 . are different from the solution of 
the linear system ([7]). Nevertheless, the resulting quadrature is expected to satisfy ([6]) with e of 
order |A„| (see Theorem IFJOl in Section B.4. 21) . since the reciprocal of i\) n can be approximated well 
by a rational function with n poles. Making the latter statement precise is the principal purpose of 
Section 0] of this paper. While the analysis of the issue is somewhat detailed, the principal idea is 
simple enough to be presented in the next few sentences. 

If P is a polynomial with m simple roots z%, . . . , z rn in (—1, 1), then the function z — > P(z)^ 1 is 
meromorphic in the complex plane; moreover, 

1 _ " 1 

p{z)~f? 1 p f (z i )-(z-z j y (15) 
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for all complex z different from zi, . . . , z m (see Theorem [57] in Section |2~5)) . The right-hand side of 
(|T5|l is referred to as "partial fractions expansion of P _1 ". Similarly, the function z — > ijjn(z)^ 1 is 
meromorphic; however, it has infinitely many poles, all of which are real and simple (see Corollary [3] 
in Section RL.l.lj) . and exactly n of which lie in (—1, 1) (see Theorem[T]m Section l2~Tj) . Suppose that 
the roots of tp n in (—1, 1) are denoted by t\ < ■ ■ ■ < t n . Motivated by (fT5|) . we analyze the partial 
fractions expansion of ^„ • ^ turns out that 

1 n 1 

*M= ^)-<t- t] > +0 ^ (16) 

for real —1 < t < 1 (see Section |4~31 and Theorem 1571 in Section I2"!!]) . In other words, (TTrj|) means 
that the reciprocal of t\j n differs from a certain rational function with n poles by a function, whose 
magnitude in the interval (—1,1) is of order |A„|. 

A rigorous version of (fTB"]) is established and proven in Section |4~31 The relation between ([5]), 
(TT2"j) . (TUT) and (|16p is studied in Section fOl The results of these two sections rely on the machinery, 
developed in Sections 14.11 14.21 



1.3 Overview of the Analysis 

1.3.1 Partial Fractions Expansion of l/ip n 

To establish (fTB]) . we proceed as follows. Suppose that x\ < X2 < • • • are the roots of ip n in (1, oo) 
(see Corollary [3] in Section [4.1 . 1[) . Suppose also that M > 1, and R > 1 is a point between xm and 
xm+i- In other words, 

1 < x\ < x 2 < ■ ■ ■ < xm < R < %m+i < (17) 

Then, for all real -1 < t < 1, 

1 " 1 

y — 



J 

^ {tf&k) ■ (t - x k ) + ^ n {-x k ) -(t + x k )) + 2^ £ R ^„(z) -{z-ty (18) 

where Tn is the boundary of the square [—R,R] x [— iR. iR] y traversed in the counterclockwise 
direction (see Theorem |2"7I in Section f^THj) . 

Suppose now that x > 1 is a root of tp n . We observe that ipn is a holomorphic function defined 
in the entire complex plane. We use the integral equation (|57|) in Section 12.11 and Theorem [55] in 
Section [2~51 to show that 



'\1> n (x + it)\* + \il/ n (x + it)\* + 1 



|C 2 • (X + it) 2 - Xn\ 

■ |</>„(1)| • %/2 



i oo (19) 



ct ■ \X n \ 

(see Theorem [3S] in Section 14.2.21) . On the other hand, we use the differential equation (1481) in 



■5 



Section |2~T1 and Theorem |2"21 in Section |2~51 to show that 



J^ n ( x + it)\^ + \^ n (x + it)\ 2 
; 1 / 4 -e ct -|^(x)|-(x 2 -l) 3 / 4 

Ct ■ (x 2 - ( X n/C 2 )) 1/4 



\(x + it) 2 -l\ 
\c 2 ■ (x + it) 2 - Xn| 



< 



(20) 



(see Theorems 1571 155 1 155 1 HP" ! 1521 in Section 11X2"]) . We combine ([TSJ and (HO]) to establish the 
inequality 



1 



< e 1 / 4 - |A„| • 



{x 2 - 1)4 



" ' (x 2 -( X „/c 2 ))t 

(see Theorem [43] in Section l4.2.2[) . Then, we use (f2Tj) to show that, for every integer M > 1, 

1 



M 

E 



< 5 



|A„|- (log(2-z M ) + (x n ) 1/4 



(21) 



(22) 



(see Theorems [44j [45] in Section [433] for a more precise statement). 

We observe that (1221) provides an upper bound on the first summand in right-hand side of (|18l) . 
While this bound is of order |A„| for xm < 0(|A„| _1 ), it diverges if we let M go to infinity (see, 
however, (|2~4"|) below). 

To overcome this obstacle, we use the integral equation (l4"4l) in Section l2~Tl to analyze the behavior 
of ip n (x) and ip' n (x) for x > |A ra | -2 (see Section f4. 3. 21) . In particular, if x > 1 is a root of ip n and if 
x > |A„|~ 2 , then 



2^n(l) 



A r , a; 



1 + b-A 



(23) 



(see Theorem [HI] in Section 14.3.21 for a more precise statement). More detailed analysis reveals that, 
if y > x > |A n |~ 2 are two consecutive roots of ip n and — 1 < i < 1 is a real number, then 



1 



V n {x)-{x-i) V n {y)-{y-t) 



' y ds 

<20 c / % 



(24) 



(see Theorem [521 in Section |4~3^]) . 

In Theorem [531 of Section 14.3.31 we establish, for all real — 1 < t < 1, the inequality of the form 



1 



fc=l 



< const • I A„ I • lo 



1 



(Xn 



,1/4 



(25) 



where (JUJ) , ([2~4"]) are used to bound the head and the tail of the infinite sum, respectively. 

Eventually, we analyze the behavior of ip n of the complex argument to demonstrate that, for all 
real — 1 < t < 1, 



lim sup 

k— ^oo 



1 



2? ™ JT Rk Tpn(z) -(z-t) 



< 2V2- |A„ 



(26) 



where is a certain sequence that tends to infinity, and the contours Tn k are as in (|18p (see 

Theorems [Ml [551 in Section 14331 for more details). We substitute (J55J and (EH) into {T5} to obtain, 



6 



for all real — 1 < t < 1, an inequality of the form 
1 A 1 



a2 m 



£jVn(*i)-(*-*j) 



< const • |A„| • ( log ( jj— 



(Xn) 



1/4 



(27) 



(see Theorems [551 EH1 m Section 14.3. 3p . In the next subsection, we overview the implications of (|2T|) 
to the analysis of the quadrature, discussed in Section [i~2l 



1.3.2 Quadrature Weights 

Roughly speaking, (l27l) asserts that, for all real —!<<<!, 



1 " 1 

M"^;fe)'(i-* 3 r 0(|Anl) 



(28) 



In other words, the left-hand side of (l28l) is uniformly bounded in (—1, 1), and its magnitude is of 
order |A„|. If we multiply both sides of (|2"8|) by ip n {t) and use (fl"2"j). we obtain 



1 = pi(t) + • • • + (p n (t) + rp n (t) ■ O (|A„|) 



(29) 



In other words, tp\ , . . . , ip n constitute a partition of unity in the interval (— 1, 1), up to an error of 
order \X n \. We integrate both sides of (j2"9"j) over (—1, 1) and use Theorem[T]in Section I2TT1 and ([T4")) 
in Section [TT21 to obtain 

2 = W 1 + --- + W n + 0(\\„\) (30) 

(see Section H. 4. 41 for more details). 

Suppose now that m ^ n is an integer. We multiply both sides of (|29p by ip m to obtain 

n 

^m(t) = M*) ■ <Pj(t) + i> m (t) ■ • O (|A„|) . (31) 

J'=l 

On the other hand, for every integer j = 1, . . . , ra, we use integration by parts to evaluate 



|A TO | 2 ■ ^ m (fj) 

-|A„| 2 



|A„ 



Wi 



ic\ n 



j) JO 



ip n {x) ■ e" 



dx 



(32) 



(see Theorem I5U1 in Section 14.4.11) . We combine (j2"7| . (|3"TT) and (I3"2l with some additional analysis to 
conclude that, for all integer < m < n, 



/l n 

1 .7 = 1 



' • |A„ ■ ( log p^— j" +Xn 



(33) 



(see Theorems |M1 Ml in Section IIX2"]) . 

According to (l3"3")l . the quadrature error © in Section H~T1 is roughly of order |A„|. It remains to 
establish for what values of n this error is smaller than the predefined accuracy parameter e > 0. In 
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Section 14.4.31 we combine Theorems 151 171 ITT1 with (|33p to achieve that goal. Namely, we show that, 
if 



2c / 1 

n > h const • log(c) • log(c) + log - 

7r \ e 



(34) 



then 



3=1 



< e, 



(35) 



for all integer < m < n (see Theorem 1(551) . 

Numerical experiments seem to indicate that the situation is even better in practice: namely, to 
achieve the desired accuracy it suffices to pick the minimal n such that |A„| < e, which occurs for 
n — 2c/tt + (9((logc) • (— loge)) (see Section [6l in particular, Conjecture [2] and Experiment 14 in 
Section [OH])- 



2 Mathematical and Numerical Preliminaries 

In this section, we introduce notation and summarize several facts to be used in the rest of the 
paper. 

2.1 Prolate Spheroidal Wave Functions 

In this subsection, we summarize several facts about the PSWFs. Unless stated otherwise, all these 
facts can be found in [3H], [3U], [H], [31], [UJ, [25], [55]. 

Given a real number c > 0, we define the operator F c : L 2 [— 1, 1] — > L 2 [— 1, 1] via the formula 

F c [if] (x) = J' <p(ty cxt dt. (36) 

Obviously, F c is compact. We denote its eigenvalues by Ao, Ai, . . . , A Tl , . . . and assume that they 
are ordered such that |A„| > |A rl +i| for all natural n > 0. We denote by ip n the eigenfunction 
corresponding to A n . In other words, the following identity holds for all integer n > and all real 
-1 < x < 1: 

A„Vn(x) = J Mty cxt dt. (37) 

We adopt the convention^ that H^nlU^-i,!.] = 1- The following theorem describes the eigenvalues 
and eigenfunctions of F c . 

Theorem 1. Suppose that c > is a real number, and that the operator F c is defined via (|36|) 
above. Then, the eigenfunctions ipo,tpi, ■ • • of F c are purely real, are orthonormal and are complete 
in L 2 [— 1 , 1] . The even-numbered functions are even, the odd-numbered ones are odd. Each function 
ip n has exactly n simple roots in (—1,1). All eigenvalues X n of F c are non-zero and simple; the 
even-numbered ones are purely real and the odd-numbered ones are purely imaginary; in particular, 

1 This convention agrees with that of 38 , 30 and differs from that of 34 . 
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We define the self-adjoint operator Q c : L 2 [— 1, 1] — > L 2 [— 1, 1] via the formula 

p -i . . 1 f 1 sin(c(x — t)) , 

Q c [<p] {x = - / ' <p{t) dt. (38) 

■K J_X X-t 

Clearly, if we denote by 2f : L 2 (R) — > L 2 (M.) the unitary Fourier transform, then 

Qc M (x) = X[-i,i](x) • J' 1 [x[-cc](0 • ?M (0] (39) 
where X[-a,a] : K — >■ K is the characteristic function of the interval [—a, a], defined via the formula 

( ) < ^ _ a — x — a ' (40) 
' a ' a ' (0 otherwise, 

for all real x. In other words, Q c represents low-passing followed by time-limiting. Q c relates to F c , 
defined via by 

Qc = f - • F* c ■ F c , (41) 
and the eigenvalues /i„ of Q n satisfy the identity 

Mn = ^-|A„| 2 , (42) 

for all integer n > 0. Obviously, 

Mn < 1, (43) 

for all integer n > 0, due to ([59]). Moreover, Q c has the same eigenfunctions as F c . In other 
words, 

, . . 1 f 1 sin (cfx — t)) , , , , .... 

for all integer n > and all — 1 < x < 1. Also, Q c is closely related to the operator P c : L 2 (U) — > 
i 2 (R), defined via the formula 

„ r i / x 1 f 00 sin (c fx — t)) , . , , , . 

P c y> x = - / — T — ¥>(* 45 

7T ./_«, X - t 

which is a widely known orthogonal projection onto the space of functions of band limit c > on 
the real line R. 

The following theorem about the eigenvalues \i n of the operator Q c , defined via (|38|) . can be 
traced back to [18] : 

Theorem 2. Suppose that c > and < a < 1 are positive real numbers, and that the operator 
Q c : L 2 [—1, 1] — > L 2 [—1, 1] is defined via (|38[) above. Suppose also that the integer N(c,a) is the 
number of the eigenvalues fx n of Q c that are greater than a. In other words, 

N(c, a) = max{fc = 1,2, . . . : [i k -i > a} . (46) 

Then, 

N(c, a) = - + ( \ log — ) logc + O flogc) . (47) 



9 



According to (|4"T1) . there are about 2c/ 'n eigenvalues whose absolute value is close to one, order of 
log c eigenvalues that decay exponentially, and the rest of them are very close to zero. 

The eigenfunctions ip n of Q c turn out to be the PSWFs, well known from classical mathematical 
physics |24j . The following theorem, proved in a more general form in [35], formalizes this statement. 

Theorem 3. For any c > 0, there exists a strictly increasing unbounded sequence of positive numbers 
Xo < Xi < • • • such that, for each integer n > 0, the differential equation 



has a solution that is continuous on [—1,1]. Moreover, all such solutions are constant multiples of 
the eigenfunction ip n of F c , defined via (|36[) above. 

Remark 1. For all real c > and all integer n > 0, (|37[) defines an analytic continuation of tp n 
onto the entire complex plane. All the roots of ip n are simple and real. In addition, the ODE (|48p 
is satisfied for all complex x. 

Many properties of the PSWF tpn depend on whether the eigenvalue Xn of the ODE (j48j) is greater 
than or less than c 2 . In the following theorem from [25], [55], we describe a simple relationship 
between c, n and Xn- 

Theorem 4. Suppose that n > 2 is a non-negative integer. 
*Ifn< (2c/ it) — 1, then x n < <? ■ 

• If n > (2c/ir), then Xn > <? ■ 

• If (2c/ V) — 1 < n < (2c/ it), then either inequality is possible. 

In the following theorem, upper and lower bounds on Xn in terms of c and n are provided. 
Theorem 5. Suppose that c > is a real number, and n > is an integer. Then, 



It turns out that, for the purposes of this paper, the inequality (|4T)]) is insufficiently sharp. More 
accurate bounds on Xn are described in the following three theorems (see [3S], [35], [37], [28] ). 

Theorem 6. Suppose that n> 2 is a positive integer, and that Xn > c 2 . Then, 



(1 - x 2 ) if>"(x) - 2x ■ ip'(x) + ( X n - cV) ip(x) = 



(48) 



n (n + 1) < Xn < n (n + 1) + c 2 . 



(49) 




(50) 



where the function E : [0, 1] — > R is defined via (|10£>[) in Section \2.S\ . 



Theorem 7. Suppose that n is a positive integer, and that 




(51) 



10 



for some 



4 

Then, 

„2 4 



0<S<^--c. (52) 



X„ > cf + - • <J • c. (53) 

7T 



Theorem 8. Suppose that n is a positive integer, and that 



2c 2c 2 p , /4 e7 rc\ 

_ < „ < _ + S ■ log — E- - 3, (54) 



7T 7T 7T 2 V 5 

for some 



4 

Then, 



3<5<^-c. (55) 



X« < c 2 + - ■ <S • c. (56) 



7T 



The following theorem is a direct consequence of Theorem [6l 

Theorem 9. Suppose that n > is a positive integer, and that 

2c , , 

n > — + 1. (57) 

X„ > c 2 + 1. (58) 



Proof. It follows from (|5T)|) of Theorem |B] that 



2c r 1 I Xn ~ c i i 



2c 2 , ? f 1 eft 2c , . . 

<— + -■ VXn - C 2 • / = — + VXn - C 2 . (59) 

7T 77 Jo Vl-i 2 'f 



We combine J59]) with <[57J) to obtain (|58jl, ■ 
In the following theorem from |27j . |28j . we provide an upper bound on |A n | in terms of n and c. 
Theorem 10. Suppose that c > is a real number, and that 

c > 22. (60) 
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Suppose also that 5 > is a real number, and that 



3 < S < — . 
16 

Suppose, in addition, that n is a positive integer, and that 

2c 2 r , f Aenc 
n> h — ■ 5 ■ log — — 

7T 7T Z \ 

Suppose furthermore that the real number £(n, c) is defined via the formula 

5 



£ (n, c) = 7056 • c • exp 



-Mi 



2nc 



Then, 



| A„ | < £(n,c). 



(61) 



(62) 



(63) 



(64) 



In the following theorem from [27] . |28] . we provide another upper bound on |A r , 
Theorem 11. Suppose that n > is a positive integer, and that 

2c rrrr 

n> h V42. 

Suppose also that the real number defined via the formula 

Xn 



Then, 



|A„| < 

1195 • c- (x n )i ■ (x n - l)i • I x r , 



exp 



7T 

4 



(65) 



(66) 



(67) 



The following theorem is a combination of certain results from 30 and |25j . 
Theorem 12. Suppose that c > is a real number, and that Xn > c 2 . Then, 

\ < V>„(1) < n +\- 



(68) 



The following theorem appears in [25], [26] . 

Theorem 13. Suppose that n > is a non-negative integer, and that x, y are two arbitrary extremum 
points of ip n in (—1,1). If \x\ < \y\, then 

\M*)\<\Mv)\- ( 69 ) 

If, in addition, Xn > c 2 ; then 

|^n(a:)|<|^„(»)|<|^„(l)|. (70) 
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The following theorem appears in [3"2"] . 
Theorem 14. For all real c > and all natural n > 1, 

max max.\ib m (t)\ <lJn. (71) 

m<n+l |t|<l ' 

In the following theorem, we provide a recurrence relation between the derivatives of ip n of 
arbitrary order (see Lemma 9.1 in [35]). 

Theorem 15. Suppose that c > is a real number, and that n > is an integer. Then, 

(1 - i 2 ) - 4*<(t) + (Xn - c 2 t 2 - 2) ip' n (t) - 2c 2 ^„(t) = (72) 

for all real t. Moreover, for all integer k > 2 and all real t, 

(1 - t 2 ) ii k+2) (t) - 2 (k + 1) *4 fc+1 > (t) + (Xn - fc (fc + 1) - c 2 t 2 ) tf> (t) 

- c 2 kt^-^ (t) - c 2 k (k - 1) ij ( n k - 2) (t) = 0. (73) 



We refer to the roots of V>n > the roots of ip' n and the turning points of the ODE (|4"5)l as " special 
points" . In the following theorem from [25] , [25] , we describe the location of some of the special 
points. 

Theorem 16 (Special points). Suppose that n > 2 is a positive integer. Suppose also that t\ < t% < 
. . . are the roots of ip n in (— 1, 1), and that si < s 2 < ■ ■ ■ are the roots oftp' n in (—1, 1). If Xn < c 2 , 
then 



/Xn /Xn 

-1 < -Z— < si < ti < s 2 < ■ ■ ■ < tn-i < s n < t n < s n+ i < < 1 (74) 

c c 

In particular, i\) n has n roots in (— 1, 1), and ip' n has n + 1 roots in (— 1, 1). On the other hand, if 
Xn > c 2 , then 



— < -1 < h < si < t 2 < ■ ■ ■ < t n -i < s„_i < t„, < 1 < y—H. (75) 
c c 

In particular, ip n has n roots in (— 1, 1), and i/4 /ias n — 1 roo£s in (— 1, 1). 

In the following theorem, proven in |25) . |26j . we describe a relation between the magnitude of 
?/> n and ij)' n in the interval (— 1, 1). 

Theorem 17. Suppose that n > is a non-negative integer, and that the functions p, q : K — > K are 
defined via f|140[) in Section \2.6l Suppose also that the functions Q,Q : (0, min {y / x«/c, l}) — > K 
are defined, respectively, via the formulae 

Q(*) = i>n(t) + -TTV • (i> n [t)) =1Pn( t ) + - ~ 2~f2 ( 76 ) 

q\l) Xn C I 



Q(t) = p(t) ■ q(t) ■ Q{t) 

= (1 - t 2 ) ■ ((Xn - c 2 t 2 ) ■ ^ 2 n {t) + (1 - t 2 ) ■ «(t)) a ) . (77) 

Then, Q is increasing in the interval (0, min {y/Xn/c, l}) , and Q is decreasing in the interval 
(0,min{ N /xw'c, l}). 
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2.2 Legendre Polynomials and PSWFs 

In this subsection, we list several well known facts about Legendre polynomials and the relationship 
between Legendre polynomials and PSWFs. All of these facts can be found, for example, in |12j . 

EE], DQ. 

The Legendre polynomials P$, Pi, P%, . . . are defined via the formulae 

Po(t) = 1, 

Pi(t) = t, (78) 

and the recurrence relation 

{k + l)P k +i(t) = {2k + l)tP k {t)-kP k -i(t), (79) 

for all k = 1,2, ... . The even-indexed Legendre polynomials are even functions, and the odd-indexed 
Legendre polynomials are odd functions. The Legendre polynomials {-Pfc}^L constitute a complete 
orthogonal system in L 2 [—1, 1]. The normalized Legendre polynomials are defined via the formula 



Pfc(t) = Pk(t) ■ Vk + I/2, (80) 
for all k = 0, 1, 2, The L 2 [—1, l]-norm of each normalized Legendre polynomial equals to one, 



i.e. 



(P fc (i)) dt=l. (81) 

Therefore, the normalized Legendre polynomials constitute an orthonormal basis for L 2 [— 1 , 1] . In 
particular, for every real c > and every integer n > 0, the prolate spheroidal wave function 
corresponding to the band limit c, can be expanded into the series 

oo oo 

m*) = E 4 n) • Ho*) - E a i n) ■ p -( x )> ( 82 ) 

for all — 1 < x < 1, where /3q™\ /3{" , . . . are defined via the formula 

P [ k n) = / t/> n (x) -Pk(x) dx, (83) 



J-i 

and ot^ , . . . are defined via the formula 

aW^-v/fcTlA (84) 
for all k = 0, 1, 2, ... . Due to the combination of Theorem [1] in Section O with (JHU), ([52]). (JHSJ), 

(^" ) ) 2 + (M" ) ) 2 + (^" ) ) 2 + ■■• = !• (85) 
The sequence /3q™^ , Pi , . . . satisfies the recurrence relation 



(86) 



A 0fi • C 


1 + A 0)2 




= Xn 


■or 


Ai,i • 


] +A h3 




= Xn 


■? { i 




Ak,k+2 ■ 




= Xn 
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for all k = 2, 3, . . . , where Ak,k> ^4fc+2,fc, Ak,k+2 are defined via the formulae 

... 2k(fc+l)-l 2 
A k .u = k(k + 1) + — „ w /, ■ c 2 



A 



k,k+2 



A k 



k+2,k 



(2fc + 3)(2fc-l) 

(ft + 2)(* + l) 



for all fc = 0,1,2, 



(2fe + 3)^(2A: + l)(2fc + 5) 
In other words, the infinite vector (^/3q U \ f3[ n \ . . .^j satisfies the identity 



(87) 



(A - Xn I). (3{ n \...y 



= 0, 



where / is the infinite identity matrix, and the non-zero entries of the infinite symmetric matrix A 
are given via ([87]) . 

The matrix A naturally splits into two infinite symmetric tridiagonal matrices, A even and A odd , 
the former consisting of the elements of A with even-indexed rows and columns, and the latter 
consisting of the elements of A with odd-indexed rows and columns. Moreover, for every pair of 
integers n, k > 0, 



0, if k + n is odd, 



(89) 



due to the combination of Theorem Q] in Section 12.11 and (|83[) . In the following theorem (that 
appears in |38j in a slightly different form) , we summarize the implications of these observations to 
the identity (|55|) . that lead to numerical algorithms for the evaluation of PSWFs. 

Theorem 18. Suppose that c > is a real number, and that the infinite tridiagonal symmetric 
matrices A even and A odd are defined, respectively, via 



A' 



and 



A 



odd 



( A),0 ^0,2 

^2,0 ^2,2 ^2,4 

^4,2 ^4,4 ^4,6 

v '■■ '■• 



A X fi 

^3,1 ^3,3 ^3,5 

^5,3 ^5,5 ^5,7 



V 



■■) 
■) 



(90) 



(91) 



where the entries Akj are defined via (|87[) . Suppose also that the unit length infinite vector /3' n ) 6 I 2 
is defined via the formula 



-An) 



where /?q"^ , f3± , . . . are defined via 




n is even, 



odd, 



If n is odd, then 



If n is even, then 

A even ^(n) = ^ . p(n) 



(92) 

(93) 
(94) 
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Remark 2. While the matrices A even and A odd are infinite, and their entries do not decay with 
increasing row or column number, the coordinates of each eigenvector ^ n > decay superexponentially 
fast (see e.g. f38l for estimates of this decay). In particular, suppose that we need to evaluate the first 
ri + l eigenvalues XQi ■ • • iXn an d the corresponding eigenvectors . . . ,^ n ' numerically. Then, we 
can replace the matrices A even , A odd in (193p . (I94p . respectively, with their N x N upper left square 
submatrices, where N is of order n, and solve the resulting symmetric tridiagonal eigenproblem by 
any standard technique (see, for example, \37\j , J^j; see also J38\j for more details about this numerical 
algorithm). The cost of this algorithm is 0{n ) operations. 

The Legendre functions of the second kind Qq, Q\, Q 2l ■ ■ . are defined via the formulae 

Qo(*) = ^log 1 + / 



2 b l-<' 

Q 1 (t) = ^l 0g l±^-1, (95) 



and the recurrence relation 

(k + 1) Qfc+i(t) = (2k + 1) tQ k {t) - fcQ fc _i(t), (96) 

for all k = 1, 2, In particular, 

^ . , 3t 2 - 1 , 1 + t 3 

QM = —^**— t -tf> 

^ , , 5t 3 -3t , 1 + t 5 2 2 

We observe that the recurrence relation (|M|) is the same as the recurrence relation (|T9")) . satisfied by 
the Legendre polynomials. It follows from (|79l) . (|M1) . that both the Legendre polynomials Pq, Pi, . . . 
and the Legendre functions of the second kind Qo, Qi, ■ • ■ satisfy another recurrence relation, namely 

t 2 P k (t) = A k _ 2 P k -2(t) + B k P k {t) + C k+2 P k+2 (t), 

t 2 Q k (t) = A k _ 2 Qk-2{t) + B k Qk(t) + C k+2 Q k+2 (t), (98) 



for all k = 2, 3, . . . , where 



(fc + l)(* 



(2/c + 3)(2fc + 5)' 



(99) 



2fc(fc + l)-l , . 

B *= (2k + 3){2k-iy (100) 

_ fc(fc-l) 

Cfe - (2Jfc-3)(2fe-l)- (101) 

In addition, for every integer fe = 0, 1, 2, . . . , the fcth Legendre polynomial P k and the fcth Legendre 
function of the second kind Q k are two independent solutions of the second order Legendre differential 
equation 

(1 - t 2 ) ■ y"{t) - 2t ■ y'{t) + k(k + 1) • y(t) = 0. (102) 
Also, for every integer k — 0, 1, . . . and all complex z such that arg (z — 1) < 71", 

Q k (z) = \ J ^ dt (103) 

(see, for example, Section 8.82 of [T2"]). 
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Remark 3. For any real number — 1 < x < 1 and integer n > 0, we can use the three-term recur- 
rences (|79p . (|96[) £o evaluate numerically Pq(x), . . . , P n {x) and Qo(x), . . . , Q n {x) with high precision, 
in 0{n) operations (see, for example, for more details). 

2.3 Elliptic Integrals 

In this subsection, we summarize several facts about elliptic integrals. These facts can be found, for 
example, in section 8.1 in [TJ], and in pQ. 

The incomplete elliptic integrals of the first and second kind are defined, respectively, by the 
formulae 

r y dt 

F(y,k)= / - (104) 
Jo Vl-& 2 sin 2 i 



E(y, k)= f \j\ - k 2 sin 2 t dt, 
Jo 



(105) 



where < y < ir/2 and < k < 1. By performing the substitution x = smt, we can write (|104p and 
(|TU5|) as 

f sia< - v) dx 
F(y ' fc) ~io ^(l-x 2 )(l-~W)' (1 ° 6) 

E(y,k)=J o ^^—^dx. (107) 
The complete elliptic integrals of the first and second kind are defined, respectively, by the formulae 

F(k) = F g, k) = f ^ M , (108) 
V2 ' Jo Vl-fc 2 sin 2 i 

/7T \ ri 2 , 

E(k)=E(-,k) = J V 1 - fc 2 sin 2 1 dt, (109) 

for all < k < 1. Moreover, 



E (VI - fc 2 ") = 1 + (- \ + log(2) - ■ fc 2 + O (fc 4 • log(fc)) . i I 1.0) 

2.4 Oscillation Properties of Second Order ODEs 

In this subsection, we state several well known facts from the general theory of second order ordinary 
differential equations (see e.g. [23 _ ). 

The following two theorems appear in Section 3.6 of [23 in a slightly different form. 

Theorem 19 (distance between roots). Suppose that h(t) is a solution of the ODE 

y"(t)+Q(t)-y(t)=0. (Ill) 
Suppose also that x < y are two consecutive roots of h(t), and that 

A 2 < Q(t) < B 2 , (112) 
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for all x < t < y. Then, 



IT TT , _ „ 

-<y-X<-. (113) 



Theorem 20. Suppose that a < b are real numbers, and that g : (a,b) —¥ K is a continuous monotone 
function. Suppose also that y(t) is a solution of the ODE 

y"(t)+g(t)-y(t) = 0, (114) 

in the interval {a,b). Suppose furthermore that 

h<t 2 <t 3 <... (us) 

are consecutive roots ofy(t). If g is non- decreasing, then 

*a - *i > <a - *2 > *4 - *3 > • • • • (116) 

If g is non-increasing, then 

*a - *i < *s - h < *4 - *3 < • • • • (117) 

The following theorem is a special case of Theorem 6.2 from Section 3.6 in |23) : 

Theorem 21. Suppose that g\,gi are continuous functions, and that, for all real t in the interval 
(a,b), the inequality gi(t) < g2(t) holds. Suppose also that the function 4>i,4>i satisfy, for all a < 
t < b, 

#'(t)+ffi(t) ■&(*)=<), 

$(*)+«»(*) •&(<)=«. (118) 
Then, fa has a root between every two consecutive roots of fa . 

Corollary 1. Suppose that the functions fa, fa are those of Theorem \21\ above. Suppose also that 

fa{to) = fa{to), fa' {t ) = <f>' 2 {to), (119) 
for some a < to < b. Then, fa has at least as many roots in (to, b) as fa. 

Proof. By Theorem 1211 we only need to show that if t\ is the minimal root of fa in (to,b), then 
there exists a root of fa in (to,t\). By contradiction, suppose that this is not the case. In addition, 
without loss of generality, suppose that fa(t),fa(t) are positive in (to,ti). Then, due to (I118p . 



and hence 



b"fa - $2 fa = (92 ~ 9i) fafa, (120) 



< / (52(a) - gi(s)) fa(s)fa(s)ds 

= [#(*)&(*)-&(*)&(*)]£ 

= 4>'i(ti)fa(ti)<0, (121) 



which is a contradiction. 
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2.5 Growth Properties of Second Order ODEs 

The following theorem appears in |19j in a more general form. We provide a proof for the sake of 
completeness. 

Theorem 22. Suppose that a < b are real numbers, and that the functions w,u,0,~f : (a,b) — > C 
are continuously differentiable. Suppose also that, for all real a < t < b, 



w'(t) 
«'(*) 



(3(t)\ (w(t) 



7(*) 



*(*) 



id that 



/3(t)^0, tWt^O, 



for all a < t < b. Suppose furthermore that the functions R, Q : (o, b) 
via the formulae 



R(t) = 



\i(t)\ 



(122) 
(123) 

are defined, respectively, 
(124) 



Q(t) = \w(t)\ 2 + R(t)-\u(t)\ 2 . 



Then, for all real a < to,t < b, 
R(t) 



R(t ) 



exp 



to 



R'(s) 

AR(s) 



|^)||7WI+»(i8(*)7(*)) 



ds 



< 



Q(t) 



Q{t 
R(t Q ) 



< 



exp 



#(*) 
4R(s) 



|/3( S )|| 7 ( S )| + 5J(/3( S ) 7 ( S )) 



ds 



(125) 



(126) 



Proof. We note that, for a each fixed t, the formula (|125p can be written in the matrix notation as 

(127) 



Q(t) = (w(t) «(*))(£ fl ( t) 



u(t) J ■ 



We differentiate Q(t) with respect to t to obtain, by using (|122[) . 

Q'(t) = w'(t)w(t) + w(t)w'{t) + R(t)u(t)u' (t) + R{t)u' (t)u(t) + R' (t)u(t)u(t) 

= /3(t)u(t)w(t) + j3(t)u(t)w(t) + R(t)~f(t)w(t)u(t) + R(t)^(t)w(t)u(t) + R 1 \t)u{t)u{t) 
(3{t) + R(t)j(t) \fw(t) 



u w\m+R(th(t) R'it) j\ U (t) 

Then, we define the functions x, y : (a, b) — > M via the formulae 

x(£) = u>(i), 



tf(t) = U{t) ■ y/Rtfj. 



(128) 



(129) 
(130) 
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We substitute (fT50"j) into (H"2"7|) . dT^5| to obtain 



Q(t) 

(x(t) y{t)) J p( t ) +R( t) 1( t) ^R^y 



(131) 



v iK*y |x(t)| 2 + ii/(t)| 2- 

To find the eigenvalues of the matrix in fj 13 1[) . we solve, for each a < t < 6, the quadratic equation 



A - 



R(t) 



A - 



y/W) J 



= o, 



(132) 



in the unknown A. Suppose that X\(t) < \2(t) are the roots of (|132p for a fixed a < t < b. We use 
dHH) to obtain 



A x (t) 

Due to (|T3Tj) . for all a < t < 6, 





(R'{t) 


2R(t) 


\2R(t) 


R'{t) 


'(R'(t) 


2R(t) 


\2R(t) 



2(\m\\i(t)\+^mh(t))) 



2(\m\h(t)\ + *WMt))) 



M*><f$ <*(*)• 



(133) 



(134) 



We substitute (|133l) into (I134[) . integrate it from to to i and exponentiate the result to obtain 

(29. ■ 

2.6 Priifer Transformations 

In this subsection, we describe the classical Priifer transformation of a second order ODE (see e.g. 
[23], [9]). Also, we describe a modification of Priifer transformation, introduced in [11] and used in 
the rest of the paper. 

Suppose that we are given the second order ODE 



dt 



(p{t)u'(t)) + q(t)u(t)=0, 



(135) 



where t varies over some interval / in which p and q are continuously diffcrentiable and have no 
roots. We define the function 9 : 1 — > K via 



P{t)u'{t) 
u(t) 



7(t) tan0(i), 



(136) 



where 7 : I — >• M. is an arbitrary positive continuously differentiable function. The function 9(t) 
satisfies, for all t in /, 



9'(t) = -^T sin 2 0(t) - 4^ cos 2 9{t) - 



Pit) 



lit) 



j(t) \ sin (20(f)) 
7W 



(137) 
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One can observe that if u'(t) = for t E I, then by (|1 36[) 

6(t) — kir, k is integer. (138) 

Similarly, if u(t) = for i € I, then 

9(t) = (k + 1/2) 7T, fc is integer. (139) 
The choice 7(f) = 1 in (| 136[) gives rise to the classical Priifer transformation (see e.g. section 4.2 in 



In [TT], the choice 7(f) = y/q(t)p(t) is suggested and shown to be more convenient numerically in 
several applications. In this paper, this choice also leads to a more convenient analytical tool than 
the classical Priifer transformation. 

Writing (gSj in the form of (fl35|) yields 

p(t) = t 2 - 1, q(t) = c 2 t 2 - X n, (140) 
for all real t > maxj^/x^/c, l}. The equation (|136l) admits the form 

^^ = V p(W)tan0(*), (141) 

which implies that 



^)=atan^|||||j + ™(t), (142) 

where m(t) is an integer determined for all t by an arbitrary choice at some t = to (the role of 7rm(i) 
in (|142j) is to enforce the continuity of 9 at the roots of ip n ). The first order ODE (|137[) admits the 
form (see HI], 0) 

0'(i) = -/(*)- sin (20(t))v(f), (143) 
where the functions /, v are defined, respectively, via the formulae 



and 



w(t) = 4 — wm — = 2 + c^^J ■ (145) 

Remark 4. In this paper, the variable t in (|14ip . (|142p . (|143p will be confined to the open ray 

(max{l, v ^fo/c},oo). (146) 
Nevertheless, a similar analysis is possible for t in the interval 

(-min{l, y/xn/c},mm.{l,y/x^/c}) . (147) 
The following theorem from \25f . \26f . summarizes such analysis for the case \ n > c 2 . 
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Theorem 23. Suppose that n > 2 is a positive integer, and that \n > c 2 . Suppose also that t\ , . . . , t n 
are the roots of ip n in (—1,1), and Si,...,s n —i are the roots of i])' n in (—1,1) (see Theorem ] 16] in 
Section Um}) . Suppose furthermore that the function 9 : [—1, 1] — > M is defined via the formula 

i(i — ~) • 7T, if ' t = ti for some 1 < i < n, 
(148) 
atan (-\/*r=S* ■ Urn) + m ® " n ' VMt)¥*o, 

where m(t) is the number of the roots of ip n in the interval (— l,t). Then, 9 has the following 
properties: 

• 9 is continuously differentiate in the interval [—1, 1]. 

• 9 satisfies, for all — 1 < t < 1, the differential equation 

9 , (t) = f(t)-v(t)-sm(26(t)), (149) 
where the functions /, v are defined, respectively, via (|144p , (|145[) in Section \2.6\ 

• for each integer < k < 2n, there is a unique solution to the equation 

6{t) = k-^, (150) 

for the unknown t in [—1.1]. More specifically, 

0(-l) = O, (151) 

9{U) = (i - • 7T, (152) 

6(s j )=j-ir, (153) 
9(l)=n-ir, (154) 

for each i = 1, . . . , n and each j ~ 1, . . . , n — 1. 

• For all real — 1 < t < 1, 

9'{t) > 0. (155) 
In other words, 9 is monotonically increasing. 

2.7 Numerical Tools 

In this subsection, we summarize several numerical techniques to be used in this paper. 
2.7.1 Newton's Method 

Newton's method solves the equation f{x) = iteratively given an initial approximation xq of the 
root x. The nth iteration is defined by 

X n = X n -i - — -. (156 

f {Xn-l) 

The convergence is quadratic provided that x is a simple root and xq is close enough to x. More 
details can be found e.g. in [7J. 
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2.7.2 The Taylor Series Method for the Solution of ODEs 

The Taylor series method for the solution of a linear second order differential equation is based on 
the Taylor formula 

u{x + h) =Y J ^f l h' + 0{h k+1 ). (157) 
4=0 J ' 

This method evaluates u(x + h) and u'(x + h) by using (|157l) and depends on the ability to compute 
u^\x) for j = 0, . . . , k. When the latter satisfy a simple recurrence relation like (|T3"]) and hence can 
be computed in 0(k) operations, this method is particularly useful. The reader is referred to [11] 
for further details. 

2.7.3 A Second Order Runge-Kutta Method 

We use the following second order Runge-Kutta Method, which can be found, for example, in [7]. 
It solves the initial value problem 

y{t ) = y , y'(t) = f(t,y) (158) 

on the interval to < t < to + L by computing 

ti+i =ti + h, 

ki+i = hf (U + i,yi + ki), 

y i+1 = Vi + (ki + k i+1 ) /2 (159) 



with i = 0, . . . , n and 



h = -, k = f(t ,yo)- (160) 
n 



Exactly n + 1 evaluations of / are required for this algorithm, which results in the total cost being 
0(n). The global truncation error is 0(h 2 ). 

2.7.4 Power and Inverse Power Methods 

The methods described in this subsection are widely known and can be found, for example, in [7]. 
Suppose that A is an n x n real symmetric matrix, whose eigenvalues satisfy 

|<ri| > H > M > ••• > K|. (i6i) 

The Power Method approximates o\ and the corresponding unit eigenvector in the following way. 



• Set i>o to be a random vector in R" such that ||t>o|| = V v o v o = !• 

• Set j = 1 and 770 = 0. 

• Compute Vj = Avj^i. 

• Set rjj = vJ_ l Vj. 

• Set Vj = Vj/\\vj\\. 

• If — 77 j x I is "sufficiently small", stop. 
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• Otherwise, set j = j + 1 and repeat the iteration. 

The output value rjj approximates o\ , and Vj approximates a unit eigenvector corresponding to o\ . 
The cost of each iteration is dominated by the cost of evaluating Avj-\. The rate of convergence 
of the algorithm is linear and equals to \a^\ / \o~i\, that is, the error after j iterations is of order 

(M/M) 3 ". 

Remark 5. A modification of the algorithm used in this paper defines rjj by 

i) ■ ( i) 

i = argmax{|i;j_i(fc)| : k = l,...,n}, rjj = — 3 ... . (162) 

The Inverse Power Method finds the eigenvalue Ofc of A and a corresponding unit eigenvector 
provided that an approximation a of is known such that 

|cr — <7fc| < max{|<7 — : j ' ^ k} . (163) 

Conceptually, the Inverse Power Method is an application of the Power Method on the matrix 
B = ( A — a I) 1 . In practice, B need not be evaluated explicitly and it suffices to be able to solve 
the linear system of equations 

(A - a I) Vj = Vj-i (164) 
for the unknown Vj on each iteration of the algorithm. 

Remark 6. If the matrix A is tridiagonal, the system f|l 64[) can be solved in 0(n) operations, for 
example, by means of Gaussian elimination or QR decomposition (see e.g J37jj, jTjj). 



2.7.5 Sturm Sequence 

The following theorem can be found, for example, in [57] (see also [5]). It provides the basis for an 
algorithm of evaluating the fcth smallest eigenvalue of a symmetric tridiagonal matrix. 



Theorem 24 (Sturm sequence). Suppose that 



/ai b 2 
b 2 a 2 b 3 



C 







• • • b n -i a„_i b n 
\0 b n a n J 



(165) 



is a symmetric tridiagonal matrix such that none of b 2 , ■ ■ ■ ,b n is zero. Then, its n eigenvalues satisfy 

(71(C) <•••< <7„(C). (166) 

Suppose also that CV- is the k x k leading principal submatrix of C , for every integer k = 1, . . . , n. 
We define the polynomials p-i,po, . . . ,p n via the formulae 



P-i{x) = 0, p (x) = 1 



(167) 



and 



p k (x) = det(C k - xl k ) , 



(168) 
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for k = 2, . . . ,n. In other words, pk is the characteristic polynomials of C'k. Then, 

p k {x) = (a k - x) pk-i(x) - b 2 k p k -2(x), (169) 

for every integer k = 1,2, ... ,n. Suppose furthermore, that, for any real number a, the integer A(o~) 
is defined to be the number of agreements of sign of consecutive elements of the sequence 

Po{o-),Pi{o-),.-.,p n (o-), (170) 

where the sign of pk{cr) is taken to be opposite to the sign of pk- i(o~) if Pk{o~) is zero. Then, the 
number of eigenvalues of C that are strictly larger than a is precisely A(a). 

Corollary 2 (Sturm bisection). The eigenvalue cxfc(C) of (|165[) can be found by means of bisection, 
each iteration of which costs 0(n) operations. 

Proof. We initialize the bisection by choosing xq < Ofc(C) < z/o- Then we set j = and iterate as 
follows. 

• Set Zj = (xj + Uj)/2. 

• If yj — Xj is small enough, stop and return Zj. 

• Compute Aj = A(zj) using (fTM)) and (flTO")) . 

• If Aj > fc, set Xj+i = Zj and j/j+i = yj. 

• If Aj < k, set Xj+\ = Xj and yj+\ = Zj. 

• Increase j by one and go to the first step. 

In the end |<Xfc(C) — Zj\ is at most yj — Xj. The cost of the algorithm is due to (|169l) and the definition 
ofA(cr). ■ 

2.8 Miscellaneous tools 

In this subsection, we list some widely know theorems of real analysis. 

The following theorem can be found in section 6.4 of [3] in a more general form. In this theorem, 
we use the following widely used notation. Suppose that g, h : (0, oo) — > C are complex- valued 
functions. The expression 

g(t) - h{t), t -> oo, (171) 

means that 

lim ^4 = 1. (172) 

t ^oo g(t) y ' 

Theorem 25 (Watson's Lemma). Suppose that b > 0, and that the function f : [0, b] — > K is twice 
continuously differentiable. Then, 

f{s)-e- st ds~ ffl, t^oo, (173) 
in the sense of (|17ip . In other words, 

t " b 



lim -— ■ / f(s) ■ e~ st ds = 1. (174) 

t^oo /(0) J Q 
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The following theorem appears, for example, in [5] in a more general form. 

Theorem 26. Suppose that xq is a real number, and u : R 2 — > R is a function of two real variables 
(t,x), defined in the shifted upper half-plane 

H Xo = {(t, x) : — oo < t < oo, xq < x < 00} . (175) 

Suppose also, that u is bounded in H Xa and is harmonic in the interior of H Xo . Suppose furthermore, 
that 

poo 

|u(i,a;o)| dt < 00. (176) 
Then, for all real t and x > x$, the value u(t, x) is given by the formula 

u(t, x) = i J ^ u(s, x Q ) ■ {t _ s) X 2 ~ X { °_ XQ)2 ds, (177) 
and, moreover, for all x > xq, 

u(t,x ) dt= u(t,x) dt. (178) 



The following theorem is a special case of the well known Cauchy's integral formula (see, for 
example, pH]!. 

Theorem 27. Suppose that D C C is an open bounded simply connected subset of the complex 
plane, and that the boundary T of D is piecewise continuously differentiable. Suppose also that the 
function g : C — > C is holomorphic in a neighborhood of D, and that none of the roots of g lies on 
r. Suppose furthermore that Z\, z%, . . . , z m € D are the roots of g in D, all of which are simple, and 
that z G D is a complex number such that g(z) 7^ 0. In other words, 

ze D\{zx,z 2 ,...,z m }, (179) 

Then, 

1 _ 1 1 r d( 

g(z) g'(z 3 ) ■ (z - Zj ) + 2m f r 9(0 • (C - z) ' ( } 

where § T denotes the contour integral over V in the counterclockwise direction. 

3 Summary 

In this section, we summarize some of the properties of prolate spheroidal wave functions (PSWFs), 
proved in the rest of the paper, mainly in Section 2) The PSWFs and the related notation were 
introduced in Section |2~T1 Throughout this section, the band limit c > is assumed to be a positive 
real number. 

In the following proposition, we describe the location of "special points" (roots of %(j n , roots of 
ifj' n , turning points of the ODE (|l8]l). in the case \n > c 2 . This proposition is proven in Theorem [29] 
and Corollary [3] in Section 14.1.11 (see also Theorem [If)] in Section 12.11) . It is illustrated in Figures [T] 
[5] (see Experiment 1 in Section l"6.1.1[) . 
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Proposition 1. Suppose that n > is a positive integer, and that \n > c 2 . Suppose also that 
X\ < x 2 < . . . are the roots of tj) n in (l,oo), and y\ < y% < . . . are the roots of ip' n in (l,oo). Then. 

1 < ^—^ <yi<x l <y 2 <x 2 < .... (181) 
c 

Also, ip n has infinitely many roots in (l,oo); all of these roots are simple. 

The following proposition summarizes the statements of Theorems 1311 1321 in Section 14.11 It is 
illustrated in Tables Q] [5J H 

Proposition 2. Suppose that n > is an integer, and that \ n > c 2 . Suppose also that x\ < Xi < . . . 
are the roots of in (l,oo). 

• For each integer k = 1,2, ... , 



7T L 1 . 7T / Xi — 1 

C 



• //, in addition, c > 1/5 and 

2c 1 



then 



• Also, 



n> - + — •(logc + log(16-e)), (183) 

7T Z7T 



7T 

^2 — xi > ~ X2 > ■ ■ ■ > Xk+i — Xk >■■■>— ■ (184) 



Xl - — > (185) 



• Moreover, 



The following proposition is an analogue of Proposition [2] in the case x« < c 2 . Its proof can be 
found in Theorem I3"3l in Section I4TT1 

Proposition 3. Suppose that n > is an integer, and that \n < c 2 - Suppose also that x± < x 2 < ■ ■ ■ 
are the roots ofipn in (Too)- Then, 

7T 

x 2 - xi < x 3 - x 2 < ■ ■ ■ < x k+1 - Xk < ■■ ■ < -■ (187) 

c 



The following inequality is proved in Theorem [35] in Section 14.2.11 and is illustrated in Tables HI 
[S] (see Experiment 5 in Section ^. 1.2j) . 
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Proposition 4. Suppose that n > is an integer, and that Xn > c 2 ■ Suppose also that x < y are 
two roots of ipn in (l>oo)- Then. 



x 2 -l 



K(*)l ■ < K(y)l < K(*)l 



X 2 - 1 C 2 !/ 2 - Xn 



^2 y 2 _ — 



(188) 



The following proposition summarizes Theorem 031 m Section 14.2.21 

Proposition 5. Suppose that n > is an integer, and that Xn > c 2 . Suppose also that x is a root 
of ifrn in (l,oo). Then, 



<e 1/4 -|AJ- 



(x 2 - l)t 



(x 2 ~ (Xn/C 2 )) 4 



(189) 



The following two estimates are proven, in a more precise form, in Theorem 1481 in Section [4.3.21 
They describe the behavior of ipn{x) for x > 1 and are meaningful only when x is large compared 
to |A„| _1 . 

Proposition 6. Suppose that n > is a non-negative integer, and that x > 1 is a real number. If 
n is even, then 



2V>n(l) 
Cx\ n 



If n is odd, then 



2^»(1) 

icx\ n 



sin(cx) + O 



cos(cx) + O 



x |A„| 4>„(l) 



x\\ n \lp n (l) 



(190) 



(191) 



The following proposition asserts that, in the interval (— 1, 1), the difference between the recip- 
rocal of ip n and a certain rational function with n poles is of order |A„|. This is an immediate 
consequence of Theorem [55] in Section 14.3.31 and the proof of Theorem [7T] in Section 14.4.41 



Proposition 7. Suppose that c > 30, and that n > is an even positive integer. Suppose also that 

(192) 



2c „ 
n > h 7. 

7T 



Suppose furthermore that —1 < t± < ■ ■ ■ < t n < 1 are the roots of ip n in (—1,1), and that the 
function I : (— 1, 1) — >■ R is defined via the formula 



for — 1 < t < 1. Then, 



= -J_ - V i 



|/(*)| < |A n | ■ ( 24 ■ log ( ) + 130 ■ (x - > 1/4 



(193) 



(194) 



for all real —1 < t < 1. 



28 



The following proposition is the principal analytical result of the paper. It is proven in Theo- 
rem IBS] in Section gX3] It is illustrated in Table [TJ and Figures |9] |TU1 |TT1 



Proposition 8. Suppose that c > is a positive real number, and that 

c > 30. 

Suppose also that e > is a positive real number, and that 

3 



exp 



(c - 20) 



< e < 1. 



Suppose furthermore that n > and < m < n are positive integers, and that 

2c / 3 , , , 1 , 1\ , fc 
n> —+ ( 10 +2 ■ 1 °g( c )+2 ■ l0g l) ■ l0g V2 



The 



/l n 
V'm(s) d,S - ^ *l J m(tj)Wj 
- 1 .7 = 1 



<£, 



(195) 



(196) 



(197) 



(198) 



where t\ < ■ • • < t n are the roots ofip n in (—1, 1), and Wi, . . . , W n are defined via (|14p in Section \1.2\ 



In Proposition [5] we address the accuracy of the quadrature, discussed in Section 11.21 More 
specifically, it asserts that to achieve the prescribed absolute accuracy e (in the sense of (jBJ), it 
suffices to take n of the order 2c/ir + O (log(c) • (log(c) — log(e))). 

The assumptions of Proposition [5] however, have a minor drawback: namely, e is assumed not 
to be "too small" , in the sense of (|196l) . In the following proposition, proven in Theorem |66j in 
Section 14.4.31 we eliminate this inconvenience. On the other hand, the resulting lower bound on n 
is considerably weaker than that of Proposition [5] 



Proposition 9. Suppose that c > is a positive real number, and that 

c > 30. 

Suppose also that e > is a positive real number, and that 

< e < 1. 

Suppose furthermore that n > and < m < n are positive integers, and that 

40* 



n- 1 



12 . . 4 , 1 
> c H • log(c) H — • log -. 

7TC ) 7T 7T £ 



Then, 



/i n 
^m(s) ds - V il> m {tj)W. 
1 i=l 



< e, 



(199) 
(200) 

(201) 
(202) 



where t\ < ■ ■ ■ < t n are the roots ofip n in (— 1, 1), and Wi, . . . , W n are defined via (|14p m S'ecizonQT^l 
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In the following proposition, we assert that the quadrature weights Wi, . . . ,W n are positive, 
provided that n is large enough. It is proven in Theorem [73] in Section 14.4.41 

Proposition 10. Suppose that c > is a positive real number, and that 

c > 30. (203) 
Suppose also that n > is a positive odd integer, and that 

n > ^ + 5 • log(c) • log (|) . (204) 

Suppose furthermore that W\, . . . ,W n are defined via (|14[) in Section \l.2\ Then, for all integer 
j = 1, . . . ,n, 

Wj > 0. (205) 



Numerical experiments seem to indicate that the assumptions (|204[) and that n be odd are 
unnecessary (see Remarks [T^J H31 in Section I4.4.4|) . 

4 Analytical Apparatus 

The purpose of this section is to provide the analytical apparatus to be used in the rest of the paper. 
4.1 Oscillation Properties of PSWFs 

In this subsection, we prove several facts about the distance between consecutive roots of PSWFs 
and find a more subtle relation between n and \n (see (|4"8l in Section I2.1[) than the inequality (|49l) . 
Throughout this subsection, c > is a positive real number and n is a non-negative integer. The 
principal results of this subsection are Theorems EU EH 

4.1.1 Elimination of the First-Order Term of the Prolate ODE 

In this subsection, we analyze the oscillation properties of ip n via transforming the ODE (|48p into 
a second-order linear ODE without the first-order term. The following theorem is the principal 
technical tool of this subsection. 

Theorem 28. Suppose that n > is a non-negative integer. Suppose also that that the functions 
^n,Qn ■ (l;Oo) — > R are defined, respectively, via the formulae 

*«(t) - M*) ■ Vt 2 -l (206) 

and 

Qn(t) = ^^ + ¥ ^, (207) 

for t > 1 . Then, 

K( t ) + Qn{t)-V n (t) = 0, (208) 

for allt > 1 . 
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Proof. We differentiate ^ n with respect to t to obtain 



I 



Then, using 



we differentiate ty' n with respect to t to obtain 

2t 



(209) 



Kit) =<(*)Vt^ : i+<W- 



<(*)V^i+<(t) 



2t 



■Mt) 



t 2 - 1 



ii> n (t) (t 2 -iy 



i 



l 



Vi 2 ^ 
-*„(t) • 



(t 2 -l) .<( t )+2t-^(t)- 

Tpn(t) ■ [Xn-C -t ) - ^— j- 



t 2 - 1 



C 2 • t 2 - Xr, 



(210) 



x t2 - 1 (i 2 -l)7 
To conclude the proof, we observe that (|208l) follows from (12 10|) . ■ 
Corollary 3. Suppose that n > is an integer. Then, i\) n has infinitely many roots in (l,oo). 
Proof. Suppose that Q n : (1, oo) — > M is defined via (|207p . Then, 

limQ„(t)=c 2 . (211) 

t— >-oo 

We conclude by combining (|21ip with (|208D of Theoreml2"51 above and Theorem [T9"l in Section l2T4l ■ 
The following theorem is a counterpart of Theorem [TB] in Section 12.11 

Theorem 29. Suppose that n > is a positive integer, and that Xn > c 2 . Suppose also that 
x\ < Xi < . . . are the roots of ip n in (l,oo), and y\ < y% < . . . are the roots of ip' n in (l,oo). Then, 



1 < 



< yi < xi < y 2 < x 2 < 



(212) 



Proof. Without loss of generality, we assume that 

M 1 ) > 0. 

We combine (|213p with the assumption that \n > c 2 and the ODE (|48p to obtain 



If, by contradiction to (|212p . 



<(i) = ^V^-^(i)>o. 



l < yi < 



(213) 



(214) 



(215) 
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then, due to 



1 Vi 

in contradiction to (12 14[) . Therefore, i\)' n is positive in the interval (l, ^/Xn/cy, in particular, 

X\ > 

c 

and 



(216) 



We combine (|217p and (|218l) to conclude that 

yXn 



<Vi< x\. 



(217) 



(218) 



(219) 



Suppose now that k is a positive integer, and y is a root of ip' n in the interval (xk,Xk+i)- Due to 



? 2 

C-V- Xr 



V>n(y)- 



(220) 



It follows from (I220p that ip' n has exactly one root between two consecutive roots of ipn- We combine 
this observation with (|219l) to obtain (|212j) . ■ 

In the following theorem, we describe several properties of the modified Priifer transformation 
(see Section l2~6l) applied to the prolate differential equation (|48|. 

Theorem 30. Suppose that n > is a positive integer, and that Xn > c 2 - Suppose also that 
Xi < x% < . . . are the roots of ip n in (l,oo), and y\ < yi < . . . are the roots of ip' n in (l,oo) (see 
Theorem ] 29\) . Suppose furthermore that the function 9 : [y^vT/c, oo) — > M is defined via the formula 



9(t) = 



ift=^, 

if t = Xi /or some i = 1, 2, 



(221) 



aton (~\/ ' xl-oW felly) + m W ' 7r ' otherwise , 
where m(t) is the number of the roots of if n in the interval (1, t). Then, 9 has the following properties: 

• 9 is continuously differ entiable in [y/Xn/c, oo). 

• satisfies, for all t > y/Xn/c, the differential equation 

9'{t) = f(t) - v{t) ■ sin(20(i)), (222) 
where the functions f,v are defined, respectively, via (|144[) . (|145l) in Section[ 
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for each integer k > — 1, there is a unique solution to the equation 



6(t) = k-~, (223) 



for the unknown t in W/Xn/ C i °°) ■ More specifically, 



(224) 

i-0-7r, (225) 
= (* - 1) • tt, (226) 

/or eac/i integer i > 1 . 

Proof. We combine (|212[) in Theorem [22] with (1221[) to conclude that # is well defined for all £ > 
-JXn/c. Obviously, is continuous, and the identities (1224[) . (|225[) . (|226[> follow immediately from 
the combination of Theorem |2T)1 and (|221j) . In addition, 8 satisfies the ODE (|222j) in (y/Xjl/c, oo) 
due to (H33), USD , (E3) in Section [2H 

Finally, to establish the uniqueness of the solution to the equation (|223p . we make the following 
observation. Due to (|221D . for any point t > s/Xn/c, the value 9{t) is an integer multiple of 7r/2 
if and only if t is either a root of ip n or a root of %j)' n . We conclude the proof by combining this 



observation with ([22411 . ([2251) and ([226]) . ■ 

The following theorem is illustrated in Table Q] (see Experiment 2 in Section 16.1.11) . 

Theorem 31. Suppose that n > is an integer, and that Xn > c 2 - Suppose also that X\ is the 
minimal root of ip n in (l,oo). Then, 

^~>h < 227 > 

Moreover, 



1 < 1 . c . ( Xl - >/*L\ (228) 



V ^1 " (X"/C 2 ) 

Proof. Suppose that y\ is the minimal root of ip' n in (l,oo). Due to Theorem l29l 



<yi< xl (229) 
c 

Moreover, due to (ggp in Theorem 1501 and l[2"T8"|) in the proof of Theorem [2^1 

sin(20(i)) > 0, (230) 

for all real yi<t<x x , where 6> is defined via ([2"2"T]) . We combine (g5DJ) with ([2"2"2]l . p25|) , (|2"2"B]) to 
obtain 



7T 

2 



0'(t) dt= (f(t) - v(t) • sin(20(i))) dt 
s/i ^ yi 



< £ f(t) dt = £ ^ C 2 - ^— ^ dt < c ■ (xi - yi). (231) 
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We combine with (ggp to obtain (j2"2"Tl) . It also follows from (j2"3"Tj) that 



which implies (JUS]). ■ 

The following theorem is a consequence of Theorems (3S] [3T] The results of the corresponding 
numerical experiments are reported in Tables [3] [3] (see Experiment 3 in Section 16.1. ID . 

Theorem 32. Suppose that n > is an integer, and that \n > c 2 . Suppose also that X\ < x% < . . . 
are the roots of ipn in (l,oo) (see Theorem \29\) . Then, 



ff- i^M-ii- ^'-^SiM^- <233) 



for each integer k = 1, 2, . . . . If, in addition, c > 1/5 and 



2 1 

n > -c+ — • (logc + log(16 • e)) , (234) 

7T Z7T 



7T 

2^2 — %1 > 2^3 — 2^2 > ' ' ' > Xk+1 ~ X]~ >■■■>— . (235) 

C 



Proof. Suppose that the functions ty n ,Q n : (l,oo) — > M. are those of Theorem l28l above. Suppose 
also that k > 1 is a positive integer. Then, due to (|207|l . 



c 2 4 ~ (X»/C 2 ) 2 t 2 - (Xn/C 2 ) 

4-1 ' t 2 - 1 



for all real Xk < t < Xk+i- We observe that ^>„ and ty n have the same roots in (1, oo) due to (|206D . 



and combine this observation with fl208[) of Theorem [28] and Theorem [19] of Section 12.41 to obtain 

(233) . 

Now we assume that c > 1/5 and that n satisfies (|234p . Also, we define the real number 5 via 
the formula 

S=l (237) 
We recall that c > 1/5 and combine (I234|) . fl23T[) and Theorem [7] in Section |2~T1 to conclude that 

o 



C 

Next, we differentiate Q n with respect to t to obtain 



> 1. (238) 



«,«) = _ -JL^ . _^ ({x . _ ^) ( , _!)_,). (239) 
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„2 \ 2 ,, „2 



We combine (I238P with Theorem |3"T1 to obtain 

( Xn -c 2 )(z?-l)-2> (^^j +7 r.^-^-2> 0, (240) 

and substitute (I240p into (|239[) to conclude that 

Q' n (t) > 0, (241) 
for alH > x\. Thus (|235l) follows from the combination of (|241j) and Theorem I2U1 in Section [2^41 ■ 

Remark 7. Extensive numerical experiments seem to indicate that, if Xn > c 2 ; i/ien (|235D always 
holds. In other words, the assumption (|234[) is unnecessary. 

The following theorem is a counterpart of Theorem [3H in the case \n < c 2 . 

Theorem 33. Suppose that n > is an integer, and that \n < c 2 . Suppose also that x\ < Xi < . . . 
are the roots oftp n in (l,oo). Then, 

7T 

X2 — Xi < X3 — X2 < ■ ■ ■ < Xk+1 — Xk < • • • < — • (242) 

c 

Proof. Suppose that the functions Q n : (1, oo) — ► R are those of Theorem above. We observe 
that Q n is monotonically decreasing in (1, oo), due to (|207p . Also, we observe that ^ n and i/j n have 
the same zeros in (1, oo), due to (|206[) . We combine these observations with (|208[) and Theorem [2^1 
in Section I2H to obtain (|232jl . ■ 

4.2 Growth Properties of PSWFs 

In this subsection, we find several bounds on \-ip n \ and \ip' n \. Throughout this subsection, c > is 
a positive real number and n is a non-negative integer. The principal result of this subsection is 
Theorem [331 

4.2.1 Transformation of the Prolate ODE into a 2x2 System 

The ODE (|48p can be transformed into a linear two-dimensional first-order system of the form 

Y'(t) = A(t)Y(t), (243) 

where the diagonal entries of A(t) vanish. The application of Theorem |2"21 in Section |2~51 to (|243[) 
yields somewhat crude but useful estimates on the magnitude of ip n and ip' n . The following theorem 
is a technical tool to be used in the rest of this subsection. This theorem is illustrated in Figures G2 
U] (see Experiment 4 in Section l6.1.2p . 

Theorem 34. Suppose that n > is a non-negative integer, and that the functions p, q : M — > M are 
defined via p40p in Section \2. 61 Suppose also that the functions Q,Q : (max { y/x^/c, l} , oo) — > K 
are defined, respectively, via the formulae 

Q(t) = €(t) + ^ ■ mt)f = €(t) + ( t2 - 1 H^ft)) 2 (244 ) 
q(t) cH z - x-a 

and 

Q(t) = p(t) ■ q(t) ■ Q{t) 

= (t 2 - 1) • ((A 2 - xn) ■ ^(t) + (t 2 1) • (^(t)) 2 ) . (245) 
Then, Q is decreasing in (max •[ y/Xn/c, 1} , oo), and Q is increasing in (max {^/Xn/c, 1} , oo). 
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Proof. We differentiate Q, denned via (|244[) , with respect to t to obtain 

<cw = 2 ■ *„(«) ■ <«) + - ^W) - («*» 



2 



2 - (!-* ) 



Xn ~ C 2 t 2 

Due to (gSJ in Section EH 



IW'CW- (246) 



^(*) = 1^2 ' - ' W*)> (247) 

for all -1 < t < 1. We substitute ([2471) into (l246| and carry out straightforward algebraic manipu- 
lations to obtain 

= (Xw _ 2 ^ 2)2 ' (X« + c 2 ~ 2c 2 t 2 ) • ^ n (t)f • (248) 
Obviously, for all real t > max { y / Xn/c, 1 } , 

Xn + c 2 - 2c 2 t 2 < 0. (249) 
We combine (|248l) with (|249p to conclude that 

Q'(t) < 0, (250) 

for all real t > max\y/Xn/ c i l}. Then, we differentiate Q, defined via (|245[) . with respect to t to 
obtain 

Q\t) = - 2t ■ (( X „ - ch 2 ) ■ v£(t) + (i - 1 2 ) ■ mt)f) 

+ (1 - t 2 ) ■ (-2c 2 t ■ 4> 2 n (t) + 2 • (xn - A 2 ) • Mt) ■ 

-2t ■ «(t)) 2 + 2-(l-t 2 )- ^ n (t) • < (t)) . (251) 

We substitute (|247[) into (|25ip and carry out straightforward algebraic manipulations to obtain 

Q'(t) = 2t ■ (2c 2 t 2 -xn- c 2 ) ■ ift(t). (252) 

We combine (|249l) with (I252[) to conclude that 

Q'{t) > 0, (253) 

for all real t > max {y/Xn/c, l}. We combine (|250l) and (1253[) to finish the proof. ■ 

Remark 8. VFe observe that the statement of Theorem \34\ is similar to that of Theorem [7j] in 
Section \2.1\ However, while in Theorem \lT\ the behavior of ip n and ip' n inside the interval (—1,1) is 
described, Theorem \34\ deals with (l,oo) instead. 

The following theorem follows directly from Theorem [M] It is illustrated in Tables 21 [5] (see 
Experiment 5 in Section [6.1. 2[) . 

Theorem 35. Suppose that n > is an integer, and that Xn > c 2 . Suppose also that x < y are two 
roots of ip n in (l,oo). Then, 



Wn{x)\ ■ < \<(V)\ < WnV)\ \ 2 l) 1 2, ■ V2 2 X i C2) - ( 254 ) 

V 1 V \Xn/c ) y i 
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Proof. Due to Theorem [ 



< x < y. 



(255) 



Due to Theorem IM1 the function Q : (s/Xn/c, oo) — > R, defined via (|244D . is monotonically decreas- 
ing. We combine this observation with (|255[) to obtain 



W n {x)\ 



x 2 - 1 



X 2 - (Xn/c 2 ) 



> 



K(»)l 



2/2-1 



(Xn/C 2 



(256) 



We rearrange (|256[) to obtain the right-hand side of (|254j) . Moreover, due to Theorem l34l the 
function Q : (y / Xn/c, oo) — > R defined via (|245p . is monotonically increasing. Therefore, 



Q(o:) = • (x 2 - 1) < K(*)l • (y 2 - 1) = VQ(») 



(257) 

which yields the left-hand side of (1254|) . ■ 
4.2.2 The Behavior of ip n in the Upper-Half Plane 

The integral equation Q37p provides the analytical continuation of ip n onto the whole complex plane. 
Moreover, the same equation describes the asymptotic behavior of ip n (x + it) for a fixed x as t grows 
to infinity (see Theorem [36] below). Comparison of these asymptotics to the estimate obtained with 
the help of Theorem [2H in Section l2~5l yields an upper bound on |f/> n (:r)| - at the roots of tjjn (see 
Theorem [41] below). The principal result of this subsection is Theorem [43] 

Theorem 36. Suppose that n > is an integer. Suppose also that x is a root of tp n in (l,oo). 
Suppose furthermore that the function Q : (0, oo) — > R is defined via the formula 



Q{t) = \if> n (x + it)\ 2 + \4>' n (x + it) | 



...a \{x + it) 2 



\c 2 (x + it) 2 - X n)V 
where i = y— 1. Then, using the asymptotic notation (|171[) of Section[ 

|^n(l)|V2 



ct \X n \ 



t — > oo, 



(258) 



(259) 



where X n is the nth eigenvalue of the integral operator (|37p . 
Proof. We use ([3"T[) in Section [2TT1 to obtain 

X n iP n (x + it) = [\ n { S )e lcsi > x+ ^ds= [ MsV csx e 



l ds 



Jo L 


r 2 - 


e ct 1 
J 


i - 


e ct 


c 


~c J 




e ct e 


-icx 


c 



e-^-^ds 



l ds 



V>„ {s/c~l)e lsx e- 8t ds. 



(260) 
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Since ip n (l) = (— 1)™ ip n (— 1), it follows from Theorem 1251 in Section \^E\ that 

Also, we differentiate (|260l) with respect to t to obtain 

A„V40 E + **) =*c / sip n (s)e lcsx e~ cts ds 



(261) 



„' „ — icx 

ie e 



2c 



(s/c - 1) ^„ (s/c - 1) e lxs e~ st ds. 



We combine (|262[) with Theorem 1251 in Section [2781 to obtain 

e ct |Vn(l)| 



|A„|t 



i — ^ oo. 



We substitute (|2"6"T]) and into (f235|> to obtain 



e ct l^(l)| 
|A„|ct 

e ct |^n(l)l 
\\Jct 



\(x + it) 2 -l\ {e ct \ip n (l)\ 



\c*{x + ity-Xn\ \ \K\t 



t — > oo, 



(262) 
(263) 

(264) 
(265) 



which implies (|259 



The rest of this subsection is dedicated to establishment of an upper bound on |^(x)| 1 at the 
roots of tp n . We start with introducing the following definition. 



Definition 1. Suppose that x > xq > 1 are real numbers. We define b c (x,xo) via the formula 

4 



b c (x 7 xo) = cxp 



x 2 - 1 



64 



E 



i 



with 



c y x 2 - xl . z -' i Si(x, Xq) + Sj(x,Xo) 



Si(x,x Q ) = x-x , 
S 2 (x,x q ) = x + x , 
S 3 (x,x ) = x-l, 
8h{x,xq) = x + 1. 



(266) 



(267) 



Next, we prove several technical theorems. 
Theorem 37. Suppose that x > xq > 1 are real numbers. Then 

where i = \f—l and, for any complex number z, we denote its real part by $R(z). 
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Proof. We fix a;o > 1, and view the integrand in (|268[) as a function of t and x. We denote this 
function by u(t, x). In other words, u(t, x) is a real- valued function of two real variables, defined via 
the formula 



u(t,x) = y 
Obviously, for fixed real x > xq 



(x + it) 2 — x\ 



{x + it) 2 - 1 



• 1 1 (a; + it) 2 - 1 



- 1. 



Next, we observe that 



5R 



(x + it) 2 — xl 
(x + it) 2 - 1 



lim u(t, x) = 0. 

\t\— ¥00 



= 1+ ^-1 



t 2 + 1- x 2 



{t 2 - (x 2 - \)) 1 +Ax 2 t 2 



and 



(X + U) 2 — Xq 



{x + it) 2 - 1 



{t 2 - {x 2 - l)Y +Ax 2 t 2 



We combine (|269l . (^711) and (g72| to conclude that for all x > x and t > 0, 



■1 < u(i,a;) < 



a;5 - 1 At 2 - Ax 2 + xl 



i (t 2 - (x 2 - 1))' +Ax 2 t 2 
Therefore, u(t, x) is a bounded function in the "shifted" upper-half plane 

H Xo = {(t,x) : x > x }. 



(269) 



(270) 



(271) 



(272) 



(273) 



(274) 



Next, again due to (|271[) and (|272[) . for all x > x and all real t satisfying the inequality t 2 > x 2 — 1, 



we have 



< u(t,x) < 



xl - 1 At 2 - Ax 2 + xl + 3 



! ( t 2 _ (3.2 _ + 4;c 2i2 
In particular, the function i — > u(t,xo) belongs to L 1 (M). In other words, 

\u(t, Xq)\ dt < oo. 



(275) 



(276) 



By carrying out tedious but straightforward calculations, one can verify that in H Xo , defined via 
(12741) . the function u(x,t) satisfies the Laplace's equation 



d 2 u , , d 2 u , 

— (t,x) + —(t,x)=0. 



(277) 



In other words, u(t, x) is a bounded harmonic function in the shifted upper-half plane H Xa . We 
apply Theorem l2"6l in Section |2~51 to conclude that, for all real t and x > xq, 



u(t, x) = 



1 



U(S,X ) 



X — Xq 



(t - s) 2 + (x - X f 



ds. 



(278) 
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and, moreover, for all x > Xq, 

fOO 

u(t,x ) dt = / u(t,x) dt. (279) 



We integrate the right-hand side of (|275l) by using the standard complex analysis residues technique 
to obtain the inequality 

, , xg-1 f°° At 2 - Ax 2 + x 2 + 3 

u(t, X) < — / — — — dt 

, ' ~ 8 J_ x (t 2 - (x - l)) 2 + 4x 2 t 2 

= if-%f- <*») 

16a; x z — 1 

We take the limit x — > oo in (|280l) and use (I279P to conclude that, for all x > Xo, 

u(t,x)dt<0. (281) 

On the other hand, due to (|27ip and (|272[) , u(t,x) is a non-negative function whenever t 2 > x 2 — 1 
and an increasing function for < t < \J x 2 — 1. Therefore, 



oo 



(t, x) > 2 • it(0, x) • \/x 2 - 1 



>-2-^-V^T=-2-J=L. (282, 
By taking the limit x — > oo in (|282D . we conclude that, for all x > xo, 

u(i,i)dt> 0. (283) 



Thus (gSHI follows from the combination of (gSD]) and (f2"53"j) . ■ 

Theorem 38. Suppose that x > xo > 1 are rea? numbers. We define the function R : R — >• R via 
i/ie formula 

R(t) = \((x + it) 2 - 1) • ((x + it) 2 - x 2 ) |" 1 . (284) 

Then, for all real t, 



R'{t) 



R(t) ^ t 2 + 5j{x,x ) 

where 5j(x,Xo) are defined via (|267p for all j = 1,2,3,4. Moreover, 



= -t-Y T? ft T7 ( 285 ) 



7o V^W/ 2 «5 l (x,x ) + l 5 J (x,xo) V ; 
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Proof. We observe that 



-m = dt logm= 2-dt l0 ^ 

Id 4 



3=1 



I d 

^^-r: l og\6 J (x,x ) + it\ 2 , 



2 dt 

3=1 

where 8%, 82, 83, 84 are defined via (|267[) . We note that, for any real number a, 



(287) 



(288) 



and thus (|285l) follows from the combination of (|287l) and (|288[) . Next, for any real numbers a, b > 0, 
t 2 dt 



o {t 2 + a 2 ) • (t 2 + b 2 ) 
i7r Res 



(z 2 + a 2 ) -(z 2 +6 2 )'' 
(1a) 2 (ib) 2 



+ Res 

7T 1 



(z 2 + a 2 ) • (z 2 + o 2 ) 



z = ib 



2ia{b 2 -a 2 ) 2ib(a 2 -b 2 )J 2 a + 6' 
and thus (|286p follows from the combination of (|285l) and (|289[) . 



(289) 



Theorem 39. Suppose that x > xq > 1 are real numbers, and that the function R : M — > K. is 
defined via (|284j) m Theorem \38l Suppose furthermore, that c, s > are reaZ numbers. Then, 



1 (x + zt) 2 - x 2 
V 2 (x + ii) 2 - 1 
cs + log 6 c (cc,cco), 

where b c (x,Xo) is defined via (|266p in Definition]!] 



1 ( (x+ity-xi \ ( R'(t) y 
' 1 (x + ii) 2 -i y v^wy 



df < 



(290) 



Proof. Suppose that the function u : M 2 — > M. is defined via (|269p in Theorem[371 Then the left-hand 
side of (12901) can be written as 



(x + it) 2 — x 2 , 



(x + it) 2 - 1 
c J dt + c u(t,x) dt+ 



1 / (s + rt) 2 -sg \ 



(x + it) 2 - 1 / \AcR(t) J 



dt 



Wi,^) + l) 2 +(^y) -V(u(t,aO + l)a)dt, 



(291) 



where the function i? 



V4ci?(t) / 

is defined via ((284)) in Theorem [38j Due to Theorem 1571 and ([275]) . 

c / u(f, a;) di < 0. (292) 
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Also, due to ([2"TT]) and (|2"72"|) in the proof of Theorem [571 for all real t > 0, 



(u(t,x) + l) 2 > (u(0,x) + 1)" = —„ X } . (293) 

x A — 1 



We combine (12931) with 02861) in Theorem El to conclude that 



2V(u(0,a:) + l) 2 Vo \^cR{t) J 32c ^ / ^ 2 a^o io 



a; 2 - 1 A 1 



64c 



where Si, 8%, S3, 84 are defined via (j267[) , and b c (x, xq) is defined via (|266p in Dcfinition[T] Thus (|290[) 
follows from the combination of {SOT}, (gg^ and (j2"M|) . ■ 

Theorem 40. Suppose that n > is an integer, and that \ n > c 2 . Suppose also that x is a root of 
ip n in (l,oo). Suppose furthermore that the function Q : K — > K is defined via the formula 



Q [t) = \i, n { x + it)\ 2 + w n { x + it)\ 2 ■ i 2 y ' 7; , . (295) 

|C \X -+■ 11) \ n 1 

Then, for all real t > 0, 
where b c is defined via (|266[) . 

Proof. We define the function <p : K — > C via the formula 

tp(t) = ip n (x + it). (297) 

Due to (g8]), 93 satisfies the ODE 

((a; + it) 2 - 1) • <p"{t) + 2i(x + it) ■ ip'(t) + (x n -c 2 (x + it) 2 ) ■ ip(t) = 0. (298) 
We define the functions w, u : K — > C via the formulae 

w(t) = tp{t), u{t) = ((x + it) 2 - 1) • <p'(t). (299) 
Due to (|298[) . the functions w,u satisfy the equation 



{x + itf - l| 



w'(t)\ ( f3(t)\ fw(t)\ 
u'(t) ~ U(t) J U(t) J ' 



(300) 



where the functions (3, 7 : R — > C are defined via the formulae 

0(t) = ((x + zt) 2 -!)" 1 , 1 {t)=c 2 {x + it) 2 -Xn- (301) 
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We combine Theorem |2"21 in Section 12751 with Theorem |3"51 above to conclude that, for all real t > 0, 



Q{t) / (R{t)Y , t , (_ Vx, 



v <w» s UroJ « J- (302) 

where b c is defined via (I266p . Q is defined via (|295|) . and the function i? : R — ► R is defined via the 
formula 

R(t) = (\( x + it) 2 -l\ ■ \{x + itf -(xn/c 2 )^ 1 . (303) 
Since ip n {x) = by assumption, it follows that 



Moreover, for all real t > 0, 

i?(t) _ (x 2 - 1) • (x 2 - ( X „/c 2 )) 



i?(0) \(x + it)*-l\-\(x + ity-(Xn/c 2 )\ 



< ^ '—±-j- 4 — VA/ " . (305) 

Thus pSBjl follows from the combination of ([502]) . ([5041 and ([505]) . ■ 

In the following theorem, we derive a lower bound on where a; is a root of ip n in (1, oo). 
It is illustrated in Tables El (see Experiment 6 in Section 16.1. 21) . 

Theorem 41 (A sharper bound on |^(x)| at roots). Suppose that n > is an integer, and that 
Xn > c 2 . Suppose also that x is a root of ipn in (l,°o)- Then, 

1 <^L. , y-»-' . b LV*i), (306) 



K(*)l " (x2-( Xn /c 2 ))3 
where b c is defined via (|266[) . 

Proof. We combine Theorem [36] with Theorem [40] and take t — > oo to conclude that 

e^|^(l)|V2 < |ft M (^-1)» ^ / ^\ 
ci|A„| cfc (3.2 _ ( x „/c 2 ))3 V c J 



which implies (|306[) . ■ 

The following theorem provides a bound on b c (x, y/Xn/ c )i defined via (|266[) in Definition Q] and 
used in Theorem |4"T1 

Theorem 42. Suppose that n > is an integer, and that Xn > c 2 . Suppose also that x is a root of 
ip n in (1, oo). Then, 

^x/^P) < e i/ 4) (308) 



C 



where b c is defined via (|266[) . 
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Proof. Obviously, b c (x, Xo), defined via (|266p . is a decreasing function of x for a fixed real number 
xq > 1. Therefore, for all real xq > 1, 



b c {x, x Q ) < b c (x 1 ,x Q ), (309) 

where x\ is the minimal root of ip n in (1, oo) (see also Theorem 129]) . We use (1267[) to conclude that 

4 



E 



i 



< 



1G 



ho- 



^ ^ ^) + Sj (x, ^) ^ 2 • 0*1 - (v^/ c )) *i - (Vx^/c) ' 

Also, due to (J12HJ) in Theorem [31] 



(310) 



- 1 2 

< - • c • I Xi 



V ^1 - (Xn/ c2 ) T 

We combine ([3"0<?1) . (j3"TUl) and (IBTTj) to conclude that 

8 





< exp 









64c 031 - (y/X^/c) 7T 



c • \Xi 



,1/4 



(311) 



(312) 



which implies (|308 



The following theorem is a direct consequence of Theorems @TJ This is the principal result 
of this subsection. 

Theorem 43 (A sharper bound on ^^(x)! at roots). Suppose that n > is an integer, and that 
Xn > c 2 . Suppose also that x is a root of ipn in (1, oo). Then, 



1 



<e 1 / 4 -|AJ- 



(x 2 - l)i 



(x 2 - (x«/c 2 )) 4 



(313) 



Proof. We combine Theorems [T2j EU [42] to obtain ([3T3]) . ■ 
4.3 Partial Fractions Expansion of l/^Vi 

In this subsection, we analyze the function l/?p n (z) of the complex variable z. This function is 
meromorphic with n simple poles inside ( — 1, 1) and infinitely many real simple poles ±xi,±X2, . . . 
outside (—1, 1) (see Theorems [Tl [TBI in Section [2TT1 and Theorem |2"91 Corollary [3] in Section l44.ip . For 
— 1 < t < 1, we use Theorem |2"T1 of Section [2T51 to construct the partial fractions expansion of l/ijj n (t) 
(see ifTS]) in Section fl . 3 . 1 D . Then, we establish that the contribution of the poles ±xi, ±X2, ... to 
this expansion is of order | A„ | . This statement is made precise in Theorems 1561 1581 which are the 
principal results of this subsection. 



4.3.1 Contribution of the Head of the Series ([Tg]) 

We use the results of Section l4~T1 and Section l4T2"l to bound the contribution of the first few summands 
of the series (fT8|) in Section 11.3.11 This is summarized in Theorem [45] below. In Theorem [44j we 
provide an upper bound on the contribution of two consecutive summands of (|18[) . Theorem [U is 
illustrated in Table [5] (see Experiment 7 in Section 16.1.31) . 
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Theorem 44 (contribution of consecutive roots) . Suppose thatn > is an integer, and thatx-n > c 2 . 
Suppose also that x < y are two consecutive roots of ip n in (l,oo). Then, 

i + i < e v4., A ,. r (z+v 2 dz 

{t-x)v n {x) + (t-y)C(y) " 1 1 L (, 2 -(Wc 2 )) 3/2 ' 

for all real t in the interval (—1, 1). 

Proof. Suppose that — 1 < i < 1 is a real number. To prove (|314[) . we distinguish between two cases. 



In the first case, 



1 



> 



1 



{x-t)W n {x)\ ~ {y-t)W n {y)[ 
We combine (|315p with Theorem l3"5l in Section l4.2.1l and Theorem l2"9l to obtain 
1 1 



(315) 



{t-x)V n (x) (t-y)ij' n (y) 
1 1 



< 



1 



1 



(x - t) K(x)| (y - t) K(y) | " (x - 1) \<P' n (x)\ (y - 1) |^(y)| 
1 



< 



1 



1 



IX 2 -(Xn/c 2 ) y 2 -l 



\4>'n( x )\ \x-l y-1 y x 2 - 1 y 2 -(xn/c 2 ) 



(316) 



We substitute pi3[) of Theorem [43] into (|316p and carry out straightforward algebraic manipulations 
to obtain 



1 



1 



(t-x)^ n (x) (t-y)i>> n (y) 

1/4. U I (- 2 - 1 ) 1 ' 



< 



I A, 



1 x 2 -(xn/c 2 ) y 2 -l 



{X 2 - {Xn/C 2 )Y \ X ~ 1 y- t \ X 2 -\ y 2 ~(Xn/c 2 )^ 

e 1 / 4 • |A„| • (x 2 - iy (x 2 - (xn/c 2 )Y (g(x) - g{y)) , 
where the function g : (y/Xn/c, oo) — > R is defined via the formula 



< 



9{z) 



z + 1 



J (z-l)(z 2 -( X n/c 2 ))- 

We differentiate (|318l) with respect to z to obtain 

y/(z - 1) (z 2 - ( Xn /c 2 )) ~2z*-2z 2 + 2z + 2{ Xn /c 2 ) 



(317) 



(318) 



> -- 



(z + l) 2 



(z-l) 2 {z 2 -{ Xn /c 2 )) 2 



(Z 2 - ( X n/C 2 )f 2 ■ Vi 2- !' 

We substitute ([gTS]) into ([3TT| to obtain 

1 1 



(319) 



< 



{t-x)V n {x) (t-y)i>> n (y) 

^.\\ n \.(x 2 -i)*(x 2 ~(x n /c'yr ■ j ye.) ,i-.< 
y ( z + i) 2 dz 



= 1/4 . 



I A, 



X ( Z 2-( Xn /c 2 )) 



3/2 • 



(320) 
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which establishes (|314[) under the assumption (|315p . If, on the other hand, 

1 1 



< 



{x-t)W n {x)\ {y-t)W n {y)Y 
then we combine (|321j) with Theorem [35] in Section 14.2.11 to obtain 
1 1 



(t-x)i>' n (x) (t-y)ipi n (y) 
1 1 



< 



1 



1 



(y - t) W n {y)\ (x - t) \V n (x)\ - {y + 1) \^' n (y)\ (x + 1) \^' n (x)\ 
1 1 f y + l x 2 — 1\ 1 y — x 



< 



We substitute (|3T^|) of Theorem g3] into p2"2l to obtain 

1 1 



(i-a:)V n (a:) (t-w)^ n (») 

(y 2 -i)* 



< 



,1/4 . 



= 1/4 . 



A,. 



I A, 



y-x 

(y 2 ~(xn/c*)) L . V-i 

v dz 



< 



(321) 



(322) 



(323) 



(z 2 - (x„/ c 2))r 

which establishes (|314|) under the assumption (132 1 1) . ■ 

The following theorem is a generalization of Theorem l4"4l 

Theorem 45. Suppose that n > is an integer, and that \n > c 2 ■ Suppose also that 1 < x\ < x-i < 
. . . are the roots of ip n in (1, oo), and that M > is an even integer. Then, for all real —1 < t < 1, 



1 



M 

E (T- x k ) ■ ^ n {x k ) 



< 4e x / 4 • |A„| • Mog(2 • x M ) + \J 1 + ^ 



(324) 



Proof. Due to Theorem l4"4l above. 



E 



M / 2 fa** 

< e 1/4 ■ |A„| • E 



(z + 1) 2 dz 



<e 1 / 4 -|A„|- / 

•/ a; 



(z + l) 2 dz 



Xl (z*-( X n/c 2 )) 



2\\3/2- 



We observe that 



z 2 dz 



(Z 2 (Xn/C 2 )f 2 

and combine (|326p with (1325[) to obtain 

1 

— ' (t — T.l. 

k 



log (z + \jz 2 - (x„/c 2 )j 



M 

E 



^ (i-a; fe ) -ip' n {x k ) 



<4e 1 / 4 -|A„|- log(2x Af ) + 



- (Xn/c 2 ) ' 



V^l ~ (Xn/c 2 ) 



(325) 



(326) 



(327) 
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It follows from the combination of Theorem [2U and Theorem [23] that 

%i 



V x i~ (Xn/c 2 ) 



1 



{Xn/c 2 ) 

X\ - {Xn/c 2 ) 



< \ 1 



^ 2c 



1 



>X:u 



2c IT V 7T 

and we substitute (I328p into (|327[) to conclude the proof. 
4.3.2 Contribution of the Tail of the Series (H]) 

In the following theorem, we establish an upper bound on \n in terms of |A r , 
Theorem 46. Suppose that n > is a positive integer, and that 

c > 30. 

Suppose also that 

1 



Then, 



lAnl < 



Xn - C < 



Proof. Suppose first that 



2c 2 10 , 
n < 1 x • — • log 



64e7r 



7T ' 7T 2 16 ~° V 10 

We combine Theorems H |8] in Section [2J] with ([329)l . (f330|) to conclude that 

X„ - c 2 < 10 • c 2 , 
provided that (|332j) holds. If, on the other hand, 

2c 2 10 , /64e7r\ 

then we combine (|334[) with Theorem [7] in Section 12.11 to obtain 

2 . 4 10 2 5 2 

Y n — C > — • • C = ■ C . 

A 7T 16 2tt 

Suppose now that the function / : (0,oo) x (l,oo) — >• R is defined via the formula 

7T • (y 2 - 1) • c 



/(c, 2 /) = 1195- 2 / lu - C -exp 



We differentiate (|336D with respect to c to obtain 

df, , f(c,y) ( 



it - (y 2 -1) -C 
4y 
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Also, we differentiate Q336I) with respect to y to obtain 



We define the real number yo via the formula 



and combine ([337]), Q338p , Q339P to conclude that 

^(c,y)<0, ^(c,y)< 0, (340) 

for all y > yo and all c > 8. Also, we defined the real number Co to be the solution of the equation 

/(c,J/o) = l, (341) 

in the unknown c > 8 (this solution is unique due to P40P ). We carry out elementary calculations 
to conclude that 

c < 30. (342) 

We combine (1339")h (|34T)j) . (j34T|) . (pH2"j) to conclude that 

/(c,2/)<l, (343) 

for all y > yo and all c > 30. Suppose now that n satisfies the inequality (I334p . We define the real 
number y n via the formula 

Vn = (344) 
V c z 

and combine J3JSJ, (|334l . (|3351) . (|3351) . (|3I2"1) . (|34U1) . (|334l with Theorem [□] in Section O to 
conclude that 

X|.|A n |</(c,y„)<l, (345) 
c z 

provided that (|334l holds. We combine (j32^|) . (j33Dl . (j333j) . (|334l . ([335]) . d33J), (|3431) to obtain 

(|331D . and thus conclude the proof. ■ 

According to Theorem l32l the distance between two large consecutive roots of ip n in (l,oo) is 
fairly close to n/c. In the following theorem, we make this observation more precise. 

Theorem 47. Suppose that n > is a positive integer, and that 

2c 

n > — + 1. (346) 

Suppose also that x,y are two consecutive roots of ip n in (l,oo), and that 

1 



< x < y. (347) 

An | 
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Suppose furthermore that 



\K\< Tq , 



and that 



The 



Xn-c < 



A„|' 



n < c ■ (y — x) < 7r + 



|A„|-x 2 



(348) 



(349) 



(350) 



Proof. Suppose that the functions $ n ,Q„ : (l,oo) — > K are those of Theorem |2"51 We combine 
Theorem E] of Section [23 (l2l£?f in the proof of Theorem Theorem HOI in Section [2H (plf)j) . ((M?)) 
and (|348[) to conclude that 



- <y-x. 
c 



(351) 



On the other hand, we combine Theorem [32] with (|3"47)) . (j348|) . (j349|l to obtain 

7T (Wc 2 )-1 



c • (y - x) <n ■ \ + 

< TT 



(Xn/C 2 ) ~ 1 

- (x™/c 2 ) 

TT 1 



< 7T- 



2 x 2 -(x„/c 2 ) 
2 

< 7T 



|A„|-a; 2 



(352) 



2-|A„| z 2 -l-l/|A n | 

Thus (|350| follows from the combination of (f35T"j) and ((352)) . I 

The following two theorems are direct consequences of the integral equation ((44)) in Section |2~T 



Theorem 48 (expansion of tp n {x)). Suppose that n > is a non-negative integer, and that x > 1 
is a real number. If n is even, then 



If n is odd, then 



tpn(x) 



cx\ n 



2^n(l) 
icx\ n 



sin(ca;) + 



1 



sin(c(a; - t))ip n (t)t 



A„V'n( 1 ) 7-i a;-* 



cos(ca;) + — 



1 



sin(c(x — t)) ip n (t)t 



iX n ip n {l)J_ 1 x-t 



dt 



(353) 



(354) 



Proof. We observe that 



11 t 

x — t x x ■ (x — t) ' 



(355) 
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for all real -1 < t < 1. We combine ((555)) with (|I2"]). (|H1) in Section to obtain 



2tt 



■■0n(x) = - 



7r 7-1 2i(a; — f) 



(if 



e 
2tt7 



7r y_ 1 2i(a; — f) 
e 



-dt 



X 



icx /•! ~ici 



2ni 



e lct Vn(t) 



27tj 7_i x(x — f) 



2ni j-i ^(^ — 



rff 
(if 



e iCT A ra ^„(-l) e-'^An^Cl) 

2irix 2irix 

1 f 1 sin (c(x — t)) , , , , 
' n -^ n (t)tdt. 



Due to Theorem [T] in Section |2~T 



and 



7r 7-i x(x — t) 



V„(l) = (-1)" • Vn(-l), 



|A„| 2 = (-l) n - A*. 

Thus {353J and ([551]) follow from the combination of ([5551) . ([5571) and (1555]) . 



(356) 



(357) 



(358) 



Theorem 49 (expansion of ^(a;)). Suppose that n > is a non-negative integer, and that x > 1 
is a real number. If n is even, then 



KM = 



2^(1) 
x\ n 



cos(cx) 

1 



sin(cx) 



1 



COS (c(x — f)) 1p n (t)t 



c\ n ip n {\) J_ 1 



CX A„Vn(l) 7-1 X — t 

1 sin (c(x - f)) -0n(f) (f 2 - 2xf) 
x [x — f) 2 



f/f 



dt 



(359) 



// n is odd, then 



2V> n (l) 



sin (CB) + £^£l + f cos(c(x-t))Mt)t dt+ 

cx \ n tp n [l) 7-1 x — t 



i f 1 sm(c(x - t))ip n (t) (f - 2xt) 



cX n tp n (l) J_± 



x (x — ty 



dt 



(360) 



Proof. The identities (|359l) . (|360[) are obtained, respectively, via straightforward differentiation of 



(|353|) . ([3541) of Theorem S5] with respect to x. ■ 

Remark 9. In the rest of this subsection, we will assume that n is even. The analysis for odd values 
of n is essentially identical, and will be omitted. 



Theorem 50. Suppose that n > is an even integer, that 

2c l 
n > hi, 

7T 



(361) 
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and that x,y are two consecutive roots of ip n in (l,oo). Suppose also that 



id that 



Suppose furthermore that 



|A„| < Tq , 



< x < y. 



Xn < 



and that the positive integer K(x) is defined via the formula 

K(x) — Round(^— ■ xj , 

where, for any real number a, Round{a) is the closest integer number to a. Then, 

2 

| sin(cir)| < 



|A„| • x' 



\cx — K(x) ■ 7r| < 



| A n | • x ' 



. cos(cx) > 1 



| At, | • x' 



and, moreover, for all real — 1 < t < 1, 



|sin(c • (y — tj) + sin(c ■ (x — t))\ < 
|cos(c • (y — t j) + cos(c • (x — t))\ < 



\Xn\-X^' 

2 

|A„|-x 2 



(362) 
(363) 

(364) 

(365) 

(366) 
(367) 
(368) 

(369) 
(370) 



Proof. We combine Theorems HI H"2l of Section |2~TI (j3"53"j) of Theorem gS] with ([3"6li]) , (|3"6"3"f to obtain 

] f 1 sin (c(xk — t)) ■ ip n (t) ■ t dt 



\sm.{cx k )\ 



A„^>„(1) J-l 

V2 



Xk - 1 



< 



\\n\(x k -l) 

2 



rdt 



a/3 | A„ | (xk - 1) ' 

which implies (|366[) . We observe that, for all real — ir/2 < s < it/ 2, 



(371) 



(372) 



and combine (|372l) with (|366l) to obtain (|367l) . The inequality (|368[) follows from the combination of 
(|3"B1)1) and (P67|) . Finally, both (131^ and (|3"7T)1) follow from the combination of (|3T)i"|) . (JM2), P>3"1) . 
and Theorem 071 ■ 
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Theorem 51. Suppose that n > is an even positive integer, and that x, y are two consecutive 
roots of ip„ in (l,oo). Suppose also that the inequalities (|361[) . Q362[) . Q363I) , (|364[) of Theorem \50\ 
hold, and that the integer K(x) is defined via (|365[) in Theorem\5Q Suppose furthermore that 

c> 1. (373) 

Then, 

= 2 ^ K(X) ^ • [1 - D(x)} (374) 
A n a: 

and 

2(-l) K{x) ib (I) 

= - 2( ij WJ • [1 - £>(*) + G(x)] , (375) 
A„y 

where the real numbers D{x) and G(x) satisfy, respectively, the inequalities 



\D(x)\ < j-^— (376) 

\A n \ ■ X 



24 



Proof. The proof is based on the identity (I359[) of Theorem |4T)1 First, we combine Theorems [TJ [T2] 
of Section O (l362|) and (|363l) to obtain 



1 sin (c(x - t)) tp n (t) (t 2 - 2xt) 



c\ n ip n (l) J_ x 

3\/2 
c-|A n |-(*-l) 2 



x (x — i) 
\-4> n {t)-t\dt < 



dt 



< 



c- |A„| • x 2 



(378) 



By the same token, 



1 



cA„Vn(l) 



1 sin {c{y-t))ib n (t) jf-2yt) 

i y(y-t) 2 



dt 



< 



c ■ |A„| • x 1 



(379) 



Also, we combine ((330)) of Theorem 07] and (1M21 . (IMB) . (1570)) of Theorem E0J to obtain, for all real 
-1 < t < 1, 

cos(c • (x — t)) cos(c • (y — t)) 



x — t y — t 

cos(c • (x — t)) + cos(c • (y — t)) cos(c(y — t)) ■ (x — y) 



x — t 
2 2 2 

|A„| • x 2 x — 1 x 2 



(y-t)- (x - t) 



< 



\\n\-X 2 



< 



(380) 
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We combine Theorems HI IT21 of Section [231 with (gSDJl to obtain 

I f 1 ( cos(c • (x — t)) cos(c-(y — t)) 



A„'-0n(l) J —I 

V2_ 8_ fl 

|An| ' X 2 



X — t 

\if> n (t)-t\dt < 



y-t 



< 



10 



|A„|-a; 2 



(381) 



We substitute (|555|> , J37DJ of Theorem HDJ ([375]) . ([3791 . dSHTJ into ([555]) of Theorem Si and use 
P73p to obtain 



-^-.(ar.^J+y.^j,)) 



< 



|cos(ca:) + cos(cy)| 



sin(c;r) sin(cy) 



4 + 4+10 



< 



24 



I A n | ' X | A n | • X 2 



(382) 



In addition, we observe that, similar to p78[) . p79[) . p80[) above, 

1 f 1 cos (c{y-t)) ■ Mt) 



X n tp n (l) J_i 
V2 



|A„|-(x-l) 



x — t 
\ip„(t)-t\dt < 



dt 



< 



|A„| • x 



(383) 



Finally, we substitute ([366]) . ([368]) . (|382l) and (1383]) into ([359]) of Theorem l49l to conclude the 
proof. ■ 

In the following theorem, we provide an upper bound on the sum of the principal parts of 1 /ip n 
at two consecutive roots of tp n in (1, oo) (see ([T8"]) in Section fl.3. It . 

Theorem 52. Suppose that n > is an even positive integer, and that x, y are two consecutive 
roots of tpn in (l,oo). Suppose also that the inequalities (|361[) . p62[) . p63l) , p64[) of Theorem \5U\ 
hold. Suppose furthermore that 



c > 1. 



Then, for all real — 1 < t < 1, 



i>' n (x)-(x-t) <P' n (y)-(y-t) 



' v ds 



(384) 
(385) 



Proof. Suppose that the integer K(x) is defined via P65I) in Theorem [BUI We combine (13741) . pTSJ), 
PTB]) . pTT) of Theorem [ST] to obtain 



(-1) K(3:) A„ 



(a:-t)(l-i?(i)) ( y -t)(l-£)(i) + G(a;)) 



1^1 . |x („ _ t) (1 _ D(i) + G(x)) - „ (i - i) (1 - D(a:))| 



jAnj 
xy 



\xyG(x) +t(y-x)(l- D{x)) - txG(x)\ < 



2 |A„G(*)| + 2|A " l(y X) = 2 |A„G(x)| + 2 |A„| • T ^. 



< 



(386) 
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where D(x),G(x) are those of Theorem IBTfl We combine Theorem Wf\ and Theorem 1501 to conclude 
that 



x z xy ■ (y — x) x n 
We substitute (|3"g7) into (|3"BT)1) and use flgggl to obtain ([555]) . 



(387) 



4.3.3 Bound on the Right-Hand Side of $TS} 

The following theorem is a consequence of Theorem 25] in Section 14.3.11 and Theorem [521 in Sec- 
tion |OH 

Theorem 53. Suppose that c > 1 is a real number, and that n > is a positive integer such that 



Suppose also that 



id that 



2c 1 
n > hi. 

7T 



|An| < Tq , 



Xn-C < 



iA„r 



(388) 
(389) 
(390) 



Suppose furthermore that 1 < X\ < x% < . . . are i/ie roots of ipn in (l,oo). Then, for all real 
-Kt<l, 



lim 



A' 



5Z ( ,// 



1 



1 



\ \"4>' n { x 2k-l) ■ (%2k-l - t) iy n (X2k) ■ (%2k - t) 

+ 20.c-|A„| 2 . 



< 



6-|A„|- (2.1og( — 



(391) 



Proof. We combine (|388p . (|389[) . (|390|) with Theorem H3 to select a positive even integer M such 
that 



1 2 



|A n |2 — | An |2- ( 392 ) 

We combine (|392[1 with Theorem [9] in Section |2~T1 and Theorem 1451 in Section f4. 3. II to obtain, for all 
real -1 < t < 1, 



M x 

E 7// iwi . 



k=l 



1p' n (xk) ■ (Xk ~ t) 



<6 - |A„| • log 



|A„| 5 



1 



(393) 



Next, we combine (|392j) with Remark [5] and Theorem [SD in Section 14.3.21 to obtain, for all real 
-1< t < 1, 



E 



1 



1 



fc=(M+2)/2 

— 2 



Tp' n {x 2 k-l) ■ (X2k-1 ~ t) 1p' n {X2k) ■ (%2k ~ t) 



< 



20 -C- 



'|A„|"= s ^ 

Thus follows from the combination of ((3"9"3")) and 



(394) 
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The rest of this subsection is devoted to the analysis of the boundary term of partial fractions 
expansion of l/?p n (see (|18l) in Section [T3H]). In the following theorem, we establish a lower bound 
on |^/) n (^)| for certain values of z. 



Theorem 54. Suppose that n > is an even positive number, and that 

w < 

Suppose also that k > is an integer number, and that 

, 8 c+1 
k> — . 

7T |A„| 

Suppose furthermore that the real number Rk is defined via the formula 

r 



7T 
C 



Then, for any real number y, 



\i/j n (R k +i- y)\ > 



^n(l) 



c • A„ 



cosh(cy) 

I^Rfc + * • y\ 



where i 



1 is the imaginary unit. Moreover, for any real number x, 



\i> n {x + i- R k )\ > 



^n(l) 



c ■ A r , 



cosh(ci?/c) 
I.T + i ■ R k \ 



(395) 



(396) 



(397) 



(398) 



(399) 



Proof. Suppose that x, y are arbitrary real numbers. We observe that 

| sin(c(a; + iy))\ 2 — | cosh(cy) • sin(cx) + i ■ cos(cx) • sinh(cy)| 2 
cosh(2cy) — cos(2cx) 



On the other hand, we combine (|396l) . (|397l) and (I401D to conclude that 

cos(2ci?t;) = cos(27rfc + tt) = — 1. 
We combine (|400|) and (|401l) to conclude that, for all real -1 < t < 1, 

| sin(c • (i? fc + % - t))| < | sin(c • (R k + iy))\ = cosh(cy). 
Next, we combine (j59"5|) . (j597|) . P%]) . Theorems [TJ in Section O to conclude that 

1 f 1 sin (c • (i? fe + % - t)) ^ n (i)t 



(400) 



(401) 



(402) 



AnV'n(l) .7-1 
cosh(cy) 2 
i4 ' |A~| 
, / s 2 cosh(cy) 
C0Sh(cy) ' UJ ' — " —4— 



Rk + iy -t 

1 , , / n , , cosh(cw) 

|^n(t) ■ t| dt < 



dt 



Rk 



< 
2 



< 



(403) 
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We combine P%]) , (|4lJ3) and ([3351) of Theorem gS] in Section gXU to obtain 

2 • ip n (l) ■ sin(c • (R k + iy)) 



\ip n {Rk + iy)\ > 



1 

1 - - 

4 



c- (R k + iy) ■ X n 
which implies (|398[) . On the other hand, due to (|400l) . 

-1 < 2 • | sin(c • (x + iR k ))\ 2 - cosh(2ci? fc ) < 1, 



for all real x. Also, due to the combination of (13951) and (1396 



16 



cosh(2ci? fc ) > exp ( -^—^ ) > 



,160 



(404) 



(405) 



(406) 



We combine (H0"5j) . flU]), (13^5] . ([3M1) . ABUT]) . Theorems [TJ [TJ in Section CO to conclude that, for all 
real x, 



1 sin (c • (a; + ii?^ — t)) ip n {t)t 



AnV'n(l) 
2 sin(c • (a; + ii?fc)) 



a; + — i 



dt 



< 



Rk 



< 



| sin(c ■ (x + 
8 



(407) 



We combine ([TO) . ([3061 . (|4U71) and (1335)) of Theorem US] in Section to obtain, for all real x, 

2-ip n (l)-sin(c-(x + iR k )) 



\ip n {x + iR k )\ > 



c • (x + iR k ) ■ A„ 



1_ i 



(408) 



which implies (|399 



In the following theorem, we use Theorem [54] to establish an upper bound on the absolute value 
of a certain contour integral. 

Theorem 55. Suppose that n > is an even positive number, and that (|395[) holds. Suppose also 
that k > is an integer number that satisfies the inequality (|396[) . and that the real number R k is 
defined via Q397D . Suppose furthermore that T k is the boundary of the square 



[-R k ,R k ] x [-i -R k ,i- R k 



(409) 



in the complex plane, traversed in the counterclockwise direction. In other words, T k admits the 
parametrization 



r k (s) 



Then, for all real — 1 < t < 1, 



' R k - iR k + 2isR k , < s < 1, 

R k + iR k - 2(s - l)R k , l<s<2, 

-R k + iR k -2i(s-2)R k , 2<s<3, 

-R k - iR k + 2(s - 3)R k , 3 < s < 4. 



(410) 



2m 



dz 



1p n {z) '(z — t) 



< 2V2 ■ \\ n \ ■ (1 + 2cR k ■ e- cR *) 



(411) 
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Proof. Suppose that — 1 < i < 1 is a real number. We combine Theorem [T2l in Section I2TT1 with 
(PS]) . (ISM)) . (15^71) . (gni) of Theorem [S2 to obtain 



r/i/ 



1 

27ri ip n (R k + iy) ■ (Rk +iy-t) 
dy 



< 



\^n{R k +iy)\ ■ \Rk + iy\ 



< 



ir ./„ cosh(cy) 



On the other hand, we combine Theorem QJ in Section O with (131?5j) . (JMS), P?T]) . of Theo- 

rem [54] to obtain 



1 

2^i 

Hit 



Rk 



dx 



ip n (x + iR k ) ■ {x + iR k - i) 



< 



dx 



< 



■■\Xn\V2 



T V-_R fc |^n(a; + i-Rfe)| • \x + iR k \ 7T • cosh(ci? fc ) 
4--y/2- |A„| -cR k - e - cRk 



dx < 



(413) 



We combine (14101) . (|412l) . (I413P with the observation that \ip n \ is symmetric about zero to obtain 

(HID. ■ 



We are now ready to prove the principal theorem of this section. It is illustrated in Table [5] and 
in Figures [SJ H] (see Experiment 8 in Section r6.1.3|) . 



Theorem 56. Suppose that c > 1, and that n > is an even positive integer. Suppose also that 

2c 

n > — + 1, (414) 

that 

|An| < ^, (415) 

and that 



Xn-c < 



A„|' 



(416) 



Suppose furthermore that —1 < t± < • • • < t n < 1 are i/ie roots of ip n in (—1,1), and £/iai £/ie 
function I : (—1, 1) — > M is defined via the formula 



= -J— - V i 



/or —1 <t< 1. Then, 

\I(t)\ < | A„ | 

where the real number I ma x is defined via the formula 

2 



/max = 24 ■ log 



13-(x™) 1/4 + 40-c-|A„| + 2v / 2. 



(417) 



(418) 



(419) 
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Proof. Suppose that 1 < x\ < x% < ■ ■ • are the roots of tf> n in (l,oo), and that k is an integer 
satisfying the inequality (I39€>|) in Theorem 1541 Suppose also that the real number Rk is defined via 
(|397p in Theorem [SH the contour Tk in the complex plane is defined via (|410|) in Theorem [55j and 
that in is the maximal root of %j) n in (f , oo); in other words, 

1< xx < ■ ■ ■ < x M < Rk < x M +i < ■■■■ (420) 

(We observe that ip n {R k ) ^ due to (j398)l in Theorem[54]) We combine (|1T7| . ([420]) and Thcorcm[27] 
of Section [2~751 to conclude that, for any real — 1 < t < 1, 



k=l 

If dz 



2-ki J Tk ip„(z) -(z-t) 

We combine the assumption that c > 1 with Theorem [9] in Section 12.11 to conclude that 



(421) 



We obtain the inequality (|418p by taking the limit fc — > oo and using (|421[) . (|422[) . Theorem [531 and 
Theorem [55] ■ 

Remark 10. T/ie conclusion of Theorem \ 56] holds for odd values ofn as well. The proof is essentially 
the same, and is based on Theorems \48\ [?M and obvious modifications of Theorems \54\ 1331 

Remark 11. Suppose that the function I : (—1, 1) — > M is defined via (]41 T[) . 7/n is even, then I is 
an even function. If n is odd, then I is an odd function. 

In the following theorem, we provide a simple condition on n that implies the inequality |A n | < 

0.1. 

Theorem 57. Suppose that c > 30, and that n > is an integer. Suppose also that 

2c 

n > — + 5. (423) 

7T 



Th 



en. 



\K\ < ^- (424) 



Proof. Suppose first that 

c> 200 • 7T. (425) 
We combine (|425p with ([12"]), (|g5j) in Section H?TJ to conclude that, in this case, 



. , , /27T • u n 2tt 1 , 

" A ""-V^ < VT < io- ( 426 ) 

On the other hand, suppose that 

30 < c < 200 • 7T. (427) 
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We observe that the interval [30, 200 • tt] is compact, and use this observation to verify numerically 
that, if (H2"Tj) holds, 

|Afloor(2c/7r+5)l < ™, (428) 

where, for a real number a, floor(a) is the largest integer less than or equal to a. We combine 
Theorem [T]in Section O (l428|) and (14261) to establish (|424|) . ■ 



In the following theorem, we summarize Theorems 1461 1561 1571 and Remark 1101 
Theorem 58. Suppose that c > is a real number and n > is an integer. Suppose also that 

c > 30, (429) 

and that 

n>— + 5. (430) 

Suppose furthermore that the function I : (—1,1) — > R is defined via the formula (|417[) in Theo- 
rem 1561 Then, 

\I{t)\ < |A n |-/ max , (431) 
where the real number I max is defined via the formula (|419[) in Theorem 1561 

Proof. We combine (|429j) . (I43Q|) with Theorem [57] to conclude that the inequality (|415j) holds. Also, 
we combine (|429p . (|415l) with Theorem [46] to conclude that the inequality (|416p holds. We combine 
these observations with Theorem 1561 and Remark [TU] to obtain (1431 1) . ■ 

4.4 PSWF-based Quadrature and its Properties 

In this subsection, we define PSWF-based quadratures of order n, find an upper bound on their 
error, and show that a prescribed absolute accuracy can be achieved by a proper choice of n. 
The principal result of this section is Theorem 1651 

Definition 2. Suppose that n > is a positive integer, and that 

-1< h <t 2 < ■■■ <t n < 1 (432) 

are the roots of ij) n the interval in (—1,1). For each integer j = l,...,n, we define the function 
ifj : (—1, 1) — > K via the formula 

V4(*i) - h) 

In addition, for each integer j = 1, . . . , n, we define the real number Wj via the formula 
We refer to the expression of the form 

n 

•/&) (435) 

j=i 

as the PSWF-based quadrature rule of order n. The points t±, . . . ,t n and the numbers W±, . . . , W n 
are referred to as the nodes and the weights of the quadrature, respectively. The purpose of (|435|) is 
to approximate the integral of a bandlimited function f over the interval [—1, 1]. 
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4.4.1 Expansion of tpj into a Prolate Series 

Suppose that n > is a positive integer. For every integer j = 1, .. .,n, we define the function 
(fj : (—1,1) — > M via (I433[) . In the following theorem, we evaluate the inner product (ipj,ip m ) 
for arbitrary m ^ n. This theorem is illustrated in Tables 111! 1121 Figure [7] (see Experiment 9 in 
Section ^J^. 

Theorem 59. Suppose that n > is a positive integer, and that m ^ n is a non-negative integer. 
Suppose also that 1 < j < n is an integer. Then, 



1 Mt) 
it — tj 



ijj m {t) dt = 



\X m \ 2 ^ m (t 3 ) r f 1 Mt) dt 



|A m | — |A n | 



i * *i 



ic\ n ^ n {l,tj) 



where tj is given via (|432[) in Definition^ and the complex-valued function \&„ : (—1,1) 
defined via the formula 



(436) 



Jo 



(437) 



Proof. We combine (13"?} with (|437p to obtain, for all real — 1 < y < 1, 



^ / & 

= -1 Jk=0 ai 

y 



ic(t — tj) 



dx dt — I e 

x=0 



i) n (t)e lcxt dtdx 



A„- / ^n(»)e "^ete = A„ • ^ n (y,tj). 
Jo 



On the other hand, 



(438) 



1 r y d 
Mt) / Tx 

t=-l Jx=0 ajL 



-icytj 



icx(t—tj) 



ic(t — tj ) 



dx dt = 1 I' Ml Ucy( t -t 3) _ A 
ic 7_i t — tj \ J 

1p n (t) dt 



e tcyt dt - 



dt 



(439) 



We combine (|438p and (|439p to obtain, for all real — 1 < x < 1, 

1 i/j„(t) dt f 1 ij) n (t) 



icX n e lcxt ^ n (x,t,) + e lcxt 
We combine §37$, (l437| and (j440|) to obtain 



i * *j 



i * *j 



Vm(i) dt= — 1p m [X) / — 

■1 I_ T j * m J x= -i Jt=-1 1 ~ l 'j 



e icxt dt dx = 



icX n r 1 0%, 
A m J—i dx 



^ m (x)e' cxt ^ n (x,t 3 ) dx+—\^J ^ m {x)e lcxt ' dx 



1 ^ n (t) dt 
l t — tj 



■(X, -tj)ty n (x,tj) dx + 1p m (tj) 



1 jj n {t) dt 
1 t _ *j 



(440) 



(441) 



We observe that tp n (—tj) — 0, and combine this observation with (|3"7)) in Section |2~T1 and (I437[) to 
obtain 



= 



^n{~tj) 

X„ 



1pn(t)( 



dt = * n (l,tj) - *„(-!, tj), 



(442) 
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and also 



AmV'mfe) - J ^ m {t)e iCtt 'dt = * m (l, -tj) - * m (-l, -tj). 

We combine ([352]) , Big) to obtain 

[*m(a: ) -*i)*r»(a;,ti)]i = _ 1 = *n(l,*j) (* m (l,-tj) - # m (-l,-t/)) 

= A m -0 m (^)*„(l,^). 
Also, we combine (1571) . Theorem Q] in Section HOI and (|437p to obtain 

* m (x, -tj)^^-(x,tj) dx = 

.1 OX 



1 



A 



^(aOe - ***' / e lc ^ / ^ m {t)e tcty dt dy dx 
m Jx=— 1 Jy=o 
1 /.l 



-— / "0m (*) / 1pn(x)e 
*m Jt=-X Jx=-\ 



f e *°&S+*)Vdy dx dt 

Jy=0 



, ll>m(t) / ^n(s) . 77 , rfa; * 

A„ f 1 lP m (t)Mt)dt (-l)" +m+1 A„ Z" 1 Vn(«) 



Am 7—1 * c (t + tj) 



-1 * *3 



tpm(t) dt. 



We combine Theorem Q] in Section [27T] with (14441) , (14451) to obtain 

ic\ n f 1 , , \ j / , w 

■(x, —tj) x M n (x,tj) dx = 



Am 7-i dx 
icX, 



An 



A^fejfutMj) + ( 1) " +mA " ^ il> m (t) dt 

*cA m 7—i i tj 



* C A„Vmfe)*„(M,) + 7^4 / P^T <Prn(t) dt. 

\A m | 7-i t — tj 

Finally, we recall that m ^ n and substitute (|446p into (14411) to obtain (|436[) . 



(443) 



(444) 



(445) 



(446) 



4.4.2 Quadrature Error 

For a positive integer n > 0, we define the PSWF-based quadrature of order n via (14321) . (|434[) 
in Definition [21 This quadrature is used to approximate the integral of an arbitrary bandlimited 
function / : (— 1, 1) — > C over the interval (—1, 1) (see ((4]) in Section fTTTI and (1435j) '). We refer to the 
difference 

f fWdt-^fM-Wj (447) 

3=1 

as the "quadrature error" (for integrating /). The following theorem, illustrated in Tables IT51 [TBT 
provides an upper bound on the absolute value of the quadrature error (for integrating ipm for 
arbitrary m < n). One of the principal goals of this paper is to investigate this error (see see © 
in Section ITTTj) . The results of additional numerical experiments, in which this quadrature is used 
for integration of certain functions, are summarized in Tables E3 [15] and Figures [HI EH EG (see 
Experiments 11, 12 in Section l"6.2.1j) . 
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Theorem 60. Suppose that n > and < m < n — 1 are integers. Suppose also that ti, . . . , t n and 
W\ , . . • , W n are, respectively, the nodes and weights of the quadrature, introduced in Definition [H 
above. Suppose furthermore that the real number P n , m is defined via the formula 



3 = 1 



(448) 



where the complex-valued function ^> n : (— 1, l) 2 — > C is that of Theorem \59\ above. Then, 



/l n 
1p m (s) ds-^ *Pm{tj)Wj 
- 1 .7 = 1 



< 



| Kn | 



oo Km 



AJ- 



\\m\ 



|Vm(0)| +C|P, 



n,m | I i 



(449) 



where \\I\\oo is the L°° -norm of the function I : (— 1, 1) 
Section \fjJ\ i.e. 



defined via (|417p in Theorem 1561 m 



|I|U =sup{|7(t)| : -!<*<!}, 



(450) 



Proof. Suppose that the function I : (—1, 1) — > R is defined via (|417[) in Theorem [551 in Section l4.3.3l 
We multiply (|4f7l) by ip n (t) ■ ip m (t) to obtain, for all real -1 < t < 1, 

iM*) = 53^m(*VjW + ll> m (t)ll> n (t)I(t), (451) 

i=i 

where, for each j = 1, . . . , n, the function <^ : (—1, 1) — > R is that of Definition [5] We combine (|37p . 
Theorem 151)1 Definition [U Theorem I5UI (|450[) , and integrate f|451[) over the interval (—1, 1) to obtain 

A TO V>m(0) = 



2 lM*i) + icXn-SlTT 1 + £ • IMU ( 452 ) 



where — l<£<lisa real number. We combine (|452l) with (14481) to obtain 

iA„r 



A m ^m(0) 



" / IA I 2 \ 

J^^ m (t j )W 3 +ic\ n P n . m +^l--L^j -Ml/Hoo. (453) 

Finally, we rearrange (|453[) to obtain (|449[) . ■ 

In the following theorem, we establish an upper bound on P n . m , defined via (|448[) above. This 
theorem is illustrated in Table [T4l and Figure \E\ (see Experiment 10 in Section IB. 1.41) . 

Theorem 61. Suppose that n,m are non-negative integers, and that < m < n. Suppose also that 
Xn > c 2 , and that the real number P rhm is defined via (|448[) in Theorem\6(A Then, 

c\P n , m \ < V32-n 2 . (454) 



G2 



Proof. Since \n > c 2 , the inequality 



f&(t)<ift{l)<n+-, 



(455) 



holds for all real — 1 < t < 1, due to Theorems [T2l [TBI in Section |2~T1 Therefore, 



Ji-i/ 



We combine (|456p with Theorem Q] in Section |2~T1 to obtain 



/' 

Jo 



ibl(t) dt> = — . 

VnK ' ~ 2 16 16 



(456) 



(457) 



We observe that 



dx 1 x 1 , x + 1 



and combine (|457l) and (|458[) to obtain 



(458) 



l-l/8n 



da; 1 1 - l/8n 1 j 2 - l/8n 

(1 - a; 2 ) 2 ~ 2 ' l-(l-l/8n) 2 + 4 ° g l/8n 



1 8n(8n-l) 1 



16n - 1 



log (16n — 1) < An + n < 5n. 



(459) 



Suppose that the functions Q(t), Q(t) : (—1, 1) — >• R are defined, respectively, via the formulae (|76l) . 
(1771) in Theorem [T7I in Section We apply Theorem [T71 with t = and < t < 1 to obtain 



Q(0)-Xn = Q(0) • p(0) • g(0) = Q(Q) 



> Q(t) = c 2 



(t 2 - 1) «(t)) 5 



(c 2 • t 2 - x») 
> c 2 ^(i) (1 - i 2 ) (X«/ C 2 - < 2 ) > C 2 V4(i) (1 - t 2 ) 
It follows from gSTJ), (|439l and (jlSOf that 



(l-t 2 ) (x„/c 2 -t 2 ) 



5n-Q(0)-^>Q(0)-^ 



l-l/8r. 



(1-x 2 )' 



> 



^ n (t) dt > 



16' 



which, in turn, implies that 



1 x 1C Xn 
< lbn ■ — — . 



Q(0) 



(460) 



(461) 



(462) 



Suppose now that j > n/2 is an integer, and tj is that of Definition [2] We combine (|462l) with 
Theorem ITTl in Section |2~T1 to obtain 



= Q(tj) > Q(0) > 



16n • Xn 



(463) 
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Due to Theorem Q3] in Section I^TTT for all integer < m < n and real — 1 < t < 1, 

\ip m (t)\ < 2rfi. (464) 
We combine Theorem[T]in Section l2~l"l with (|437D of TheoremlSTflabove to obtain, for all real < t < 1, 



l**(M)| 



4>„(x)e~ 



dt 



<\J_^n{x)\dx<^. (465) 



Finally, we combine (14151) . (gM|), (llMl) and dlllS]) to obtain 



c\P n , m \ < cn ■ max 



*n(l,*i) 



\/l6^ \/2 „ _ 

< cn • • — — • zyn, 

c 2 



(466) 



which implies (|454[) . 

Corollary 4. Suppose that m is an odd integer. Then, P n 



0. 



Proof. Suppose that 1 < j < n is an integer, and t\,...,t n are the roots of tp n in (—1,1). We 
combine Theorem [1] and (j3"7| in Section [2TT1 with (|437l) to obtain, for every j = 1, ...,n, 

(-l) B .* n j(l) + * n , B+ l_ i (l) = 
1 ,1 

^„(-x)e- lc ^ da;+ / i, n (x)e ,cxt ' dx = 



(467) 



We observe that V4 is odd for even n and even for odd n, and combine this observation with (|467p 
to obtain, for every integer j = 1, . . . , n, 



*n,j(l) _ *„, n+ i_j(l) 



We combine (14681) with (1448ft to obtain 



(468) 



(469) 



In the following theorem, we simplify the inequality (|449[) of Theorem [501 It is illustrated in 
Table [TBI and in Figure[9](see Experiment 12 in Section l6.2.1l) . See also Conjecture [2] and Remark l26l 
in Section ^. 2. II 

Theorem 62. Suppose that n > and < m < n ~ 1 are integers. Suppose also that ti, . . . , t n and 
W\ , . . • , W n are, respectively, the nodes and weights of the quadrature, introduced in Definition 
above. Suppose furthermore that 



c > 30, 



(470) 
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id that 



2c r 
n > h 5. 

7T 



The 



/l n 
-0m (s) OS - y^V> m (tj)Wj 
1 J=l 



< |A n |- (24 -log) ^ ) + 6- Xl 



(471) 



(472) 



Proof. We combine Theorems [TJ HU in Section 12.11 the inequality (|471l) and Theorems [501 [HI] to 
conclude that 



/l n 
■0 m (s) ds - ^2ip m (tj)Wj 

- 1 .7 = 1 



< ||/|U + |A„|-(2^+V32-n 2 ) 



(473) 



where |!/||oo is dehned via (I450P in Theorem IBTil Next, we combine (|470l) . ([4711) . Theorem [9J in 
Section HTJ Theorems [571 [551 in Section |4"X31 and (jlBD]) to conclude that 



||I||oo < |A„| • (^24 • log (^-A + 13 • ( X „) 1/4 + 4^ + 2y/2j 
We combine ([4711) with Theorem [S] in Section l!01 to conclude that 



(474) 



(475) 



Also, we observe that, due to the combination of (I470P and Theorem [§J in Section [2~T1 
V32 • X n + + 15 • (x„) 1/4 + 2V2 + 24 • log(2) = 

Xn • (V32 + 4 • Xn 1/2 + 15 ' Xn 3/4 + (2^2 + 24 • log(2)) ■ x" 1 ) < 6 • X n- (476) 

Now (|4"T2"]) follows from the combination of (|4T3"1) . (I4T41 . (|4T5|) and (fTT5|) . ■ 

The following theorem is a conclusion of Theorem [TT] of Section [2TT1 and Theorems IB"2l [571 above. 

Theorem 63. Suppose that n > and < m < n — 1 are integers. Suppose also that t\,...,t n and 
W\ , . . . , W n are, respectively, the nodes and weights of the quadrature, introduced in Definition 
above. Suppose furthermore that 



and that 



Then, 



c > 30, 



2c - 
n > h 7. 

7T 



/I ri 

- 1 .7=1 



< 14340 • • exp 

c 7 



I" Xr, 



(477) 



(478) 



(479) 
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Proof. We combine (14771) . (|478j) with Theorem 1521 above to obtain 



/l n 

- 1 3=1 



;|A„|.(24-log( — )+6-Xr 



Suppose first that 



Then, 



I A n | < exp 



Xn 

4 



AJ ■ 24-loe 



1 



6-Xn <48-|An|-log 



1 



We combine (|477l) , (|481l) and Theorem 0] in Section to conclude that 



|A n | < exp 



c 



(480) 



(481) 



(482) 



(483) 



We combine (|38"T]) . (|4"g2l and (jl55|) to obtain 

1 



A,!' ( 24 • In;; ( — j + 6 • Xn j < 48 • |A„| • log ^— I 



48 • exp 

Suppose, on the other hand, that 



Xn , n 

— = 12 • Xn ■ exp 



Xn 

' 4 



exp 



Xn 

T 



< | A„ | < 



10 



(484) 



(485) 



(note that the right-hand side inequality in (|485[) follows from the combination of (|477l) . (|478[) and 
Theorem l57|) . It follows from (|485[) that, in this case, 



(486) 



|A„| ■ ( 24- log ( j— ) +6- X n ) < 12-Xn-IAr, 



We combine (|478l) with Theorem QTJ to obtain 

|A„| < 1195 - ^ - exp 



K Xn ~ (T 

4 ' 



'Xr, 



(487) 



We combine (|477l) with Theorem @] of Section 12.11 to conclude that 



exp 



Xn 

T 



X. n 



< 1195 • ^ • exp 



K Xn ~ C 

4 ' 



'Xn 



We combine (BED, (SHI, (053, 0SZD, (SEE} that 

1 



An| ' 24 -log 



|A„| 



12- x« • H95- % - exp 



6 • Xn J < 

_7T X„ - C 2 
4 ' 



(489) 



Now (14791) follows from the combination of (ggQJI and (|489j) . 
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4.4.3 The Principal Result 

In Theorem 1631 we established an upper bound on the quadrature error for integrating i\) m (see 
(|479D ). However, this bound depends on Xn- In particular, it is not obvious how large n should be 
to make sure that the quadrature error does not exceed given e > 0. In this subsection, we eliminate 
this inconvenience. 

The following theorem is illustrated in Table \W\ (see Experiment 14 in Scction r6.2.1[) . 



Theorem 64. Suppose that c > is a positive real number, and that 

c > 30. 

Suppose also that e > is a positive real number, and that 

< log - < -^jL • c - 3 • log(c) - log(6 5 • 14340). 
Suppose furthermore that the real number a is defined via the formula 



4^ 



log - + 3 • log(c) + log(6° • 14340) ) , 



and that the real number v(a) is defined via the formula 

, 2c a 
v [a) = h — • log 

7T Z7T 



16ec 



Suppose, in addition, that n > and < m < n — 1 are integers, and that 

n > u{a). 

Then, 



/i n 
4> m {s) ds-^ lp m {tj)Wj 
1 i=l 



< e, 



where tj, Wj are defined, respectively, via (|432j) . (I434[) in Definition^ 
Proof. It follows from (|4UTj) that 



5c > a > 



(3 • log(c) + log(6 5 • 14340)) , 



where a is defined via (|492p . We observe that 

/ 16ec 



log 



16c 
a 



(490) 



(491) 



(492) 



(493) 



(494) 



(495) 



(496) 



(497) 



d 

— a ■ log . 

da i \ a 

and hence the function v : (0, 16c) — > M, defined via ([4931) . is monotonically increasing. We combine 
(|4"90l) . (|4"92l) . (I496| . (|497|) to conclude that 



2c „„ , s 2c 5c , 
h 30 < via) < h — • log 

7T 7T 2ir 



16e 



< 



5c 



(498) 
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We combine Theorem [7] of Section [O with (RM)) . (|4lH|) and (|4lMl) to obtain the inequality 

Xn > c 2 + a ■ c. (499) 
Suppose now that the function / : (c, oo) — > M is defined via the formula 

/(y) = y W ■ exp 



2 ? 

7T y z — c 



4 y 

We differentiate (|500l) with respect to y and use (|490l) to obtain 



(500) 



f'iv) 



f(y) 



y 



io -y i + 



< o, 



(501) 



for all y > c. We combine (|490|) . (j498| . (|499| . (l500| . (I50H with Theorem [63] to conclude that 



/l n 

- 1 3 = 1 



X 5 

14340 • -f ■ exp 
c' 



7[ Xn ~ C 

4 ' 



< 



< 



14340 • c 3 • (l + -J 
We combine gSBJ, (j502l to obtain 

/l n 
-0m (s) (fe - 1p m (tj)Wj 
1 .7 = 1 



• exp 



4 ^1 + a/c 



< 14340 • 6 5 • c 3 • exp 



7T a 

4 ' 7i. 



(502) 



(503) 



Now (gg5| follows from the combination of (jl9"2"f and (|5U3"1) . 



The following theorem is a direct consequence of Theorem 1641 This theorem is one of the 
principal results of the paper. It is illustrated in Table [T9l fsee Experiment 14 in Section f6. 2. ip . See 
also Conjecture [5] in Section ^. 2. II 



Theorem 65. Suppose that c > is a positive real number, and that 

c > 30. 

Suppose also that e > is a positive real number, and that 

3 



exp 



(c - 20) 



< e < 1. 



Suppose furthermore that n > and < m < n are positive integers, and that 

2c / 3 , . , 1 , 1\ , /c\ 
n>-+ 10+-.log(c) + 5 .log- -logy. 



Then. 



/i n 
m (s) ds - Ipmit^Wj 
- 1 .7 = 1 



< e, 



(504) 



(505) 



(506) 



(507) 



where tj, Wj are defined, respectively, via (|432p . (14341) in Definition^ 
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Proof. We observe that, for all real x > 30, 



| • (a; - 20) < — ^ • x - 3 • log(x) - log(6 5 • 14340). 
2 4y6 



Also, we combine f|504[) . (|505[) to conclude that 

^ • (log - + 3 ■ log(c) + log(6 5 • 14340)") < 10 + | • log(c 



1 , 1 
- ■ log - . 

2 & e 



Furthermore, we combine (|504l) . (|505l) to conclude that 

— • flog - + 3 ■ log(c) + log(6 5 • 14340)) > 89 > 2 • 16e. 



(508) 



(509) 



(510) 



Now (|5U71) follows from the combination of ([504]) . ([5D5]) . pMl) . fSDS]) . (jSTU)) and Theorem |M1 



The assumptions of Theorem [53] contain a minor inconvenience - namely, the parameter s is 
not allowed to be "too small" (in the sense of (|505|) ^ . In the following theorem, we eliminate this 
restriction. On the other hand, for the values of e in the range (|505p . the resulting inequality for n 



(511) 
(512) 



is much weaker than (|506|) . 

Theorem 66. Suppose that c > is a positive real number, and that 

c > 30. 

Suppose also that e > is a positive real number, and that 

< s < 1. 

Suppose furthermore that n > and < m < n are positive integers, and that 

i 40 \ 12 ... 4 . 1 
n ■ ( 1 >cH • log(c) H — • log -. 

7TC / 7T 7T E 



The 



/l n 
VVn(s) - ^ lp m (tj)Wj 
-1 .7 = 1 



< e, 



(513) 



(514) 



where tj, Wj are defined, respectively, via (|432[) . (I434[) m Definition^ 

Proof. We combine f|513[) with Theorem [5J in Section 12.11 and (|109[) in Section 12.31 to conclude that 



c 2 <n 2 < Xr, 



(515) 



Also, we combine (jBTTj) . ([5T2"]) . (j5T3")) , (j5"T5l) with Theorem 1531 and PJT|) in the proof of Theorem [ 
to conclude that 



/l n 
%b m (s) ds - y2^Pm(tj)Wj 
- 1 3=1 



14340 • ^ • exp 

c' 



7T Xn ~ C 

4 ' 



'Xr, 



< 



< 



14340 • c d 



exp 



4 \c nl 



(516) 
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We take the logarithm of both sides of (|516[) and use (|515[) to obtain 



log 



/l n 
1 i=l 



< 



7T 

n + — • c < 
4 



-) - 1 

(Tt\ 7T 7T 
- -10 n + - ■ c < 
cJ 4 4 

/12 , , N / 40 \ \ . . 

log(c) - n ■ 1 + c . (517) 



4 \ 7T \ 7TC_ 

Now (15 14[) follows from the combination of (|513l) and (|5 1 7[) . ■ 
4.4.4 Quadrature Weights 

In this subsection, we analyze the weights W\, . . . , W n of the quadrature, defined in Definition [2] 
in Section 14.41 This analysis has two principal purposes. On the one hand, it provides the basis 
for a fast algorithm for the evaluation of the weights. On the other hand, it provides a theoretical 
explanation of some empirically observed properties of the weights. 

The results of this subsection are illustrated in Table l20l and in Figure [12] (see Experiment 15 in 
Section [5X2]). 

In the following theorem, we describe a function, whose values at the roots ti,...,t n of ip n in 
(—1,1) are equal to the quadrature weights W\, . . . , W n , up to a certain scaling. 

Theorem 67. Suppose that n is a non-negative integer. Suppose also that the function <!>„ : 
(— 1, 1) — >• K is defined via the formula 



$ n (t) = J2a£ ) Q k (t), (518) 



."'V 

fc=0 

where Qk(t) is the kth Legendre function of the second kind, defined in Section W^. and is the 
kth coefficient of the Legendre expansion ofip n , defined via (|84[) in Section \2."A Suppose furthermore 
that t\ < ■ ■ ■ < t n are the roots ofip n in (—1,1). Then, for every integer j = 1, 2, . . . , n, 

~^) = \fjf^. (519) 



Proof. Suppose that 1 < j < n is an integer, and that S > is a positive real number. We combine 
USTSJ with (P ]l . (|gg ]) . ([51 ) . (fTU31 in Section O to obtain 

provided that S is sufficiently small. Suppose now that e > is a real number, and that 

e< i-mm{fe-l|, fe + l|}. (521) 
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We observe that, since tj is a root of ipn, the right-hand side of (|519[) is well defined. We combine 
this observation with (I520p , (I521[) to evaluate 

s^o m 5>o \2 J t ._ e tj + iS-t 2J tj _ e tj-t J 

lim \ I ip n (tj + s) ■ ( — 5— ] ds = 

s^o, s>o 2 J_, J \s + iS si 



lim 



iS ■ tp' n (tj) f £ ds 



lim S ■ ip' n (tj) ■ arctan ( ~") = 0. (522) 



8-tO, 8>0 2 J_ e S + iS 

lim 8 ■ ib'(tj) ■ arctan ( — 
We combine ([515]) , ([520]) , (|522]l to obtain ([515)1 ■ ■ 

The following corollary is a direct consequence of Definition [2] and Theorem l67l 

Corollary 5. Suppose that n > and 1 < j < n are positive integers. Suppose also that the function 
<&„ : (— 1, 1) — > R is defined via (I518[) m Theorem\67\ Then, 

w >- 2 'Wty (523) 

where tj , Wj are defined, respectively, via (|432p , (I434[) o/ Definition [D 

Corollary [5] is illustrated in Table [20) We observe that Theorem [57] and Corollary [5J describe a 
connection between the weights Wi, . . . , W„ and the values of $ n at t\, . . . , t n , where the function 
$„ is defined via f)518[) . In the following theorem, we prove that $ n satisfies a certain second-order 
non-homogeneous ODE, closely related to the prolate ODE (|48l) in Section |2~T1 



Theorem 68. Suppose that n is a non-negative integer, and that the function $„ : (—1,1) — > R 
is defined via (15 1 8|) in Theorem \6l\ Suppose also that the second-order differential operator L n is 
defined via the formula 

L n [<p] (t) = (1 - t 2 ) <p"(t) - 2t<p'(t) + ( X n - c 2 t 2 ) <p(t). (524) 
Then, in the interval (—1,1) the function satisfies the nonhomogeneous ODE 



(t) = -c 2 (4 n) t + al n) /3), (525) 



where the coefficients o^\cti are the first two coefficients of the Legendre expansion ofip n , defined 



via (|84J) in Section [Ol 
Proof. We combine ([TU2"]) . ([Ml) of Section O with ([524]) to obtain 

in [Qk] = (Xn - k(k + 1) - c 2 t 2 ) ■ Q k , (526) 
where Qk is the fcth Legendre function of the second kind, defined in Section 12.21 We combine 
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of Section with fl52g} to obtain 



E 



,fe=0 



= E «i n) - + 1) - c2t2 ) Qk = 

oo 

(Xn-*(* + l))0* 

fc=0 

oo 

" c 2 E a i" } (^-2^-2 + S fc Q fc + C k+2 Q k+2 ) = 

k=0 

oo 

£ [( X „ - + 1)) a<"> - c 2 ( Q W A k + 4 l) B k + a{%C k 

+ [(Xn 1(1 + 1)) a[ n) ~ c 2 (4"^! + a^B,)] Q 



k=2 



Qo 



( Xn - 0(0 + 1)) ag> - c 2 (o4 n) A + 4 n) B 

4" } (< 2 g! - B 1 Q 1 - C 3 Q 3 ) + 4 K) (* 2 <3o - P0Q0 - C 2 Q 2 ) 



(527) 

where B k ,C k are defined, respectively, via (j99j), fllOOf) . (|101[) in Section [2T2l By the same token, 
(|527l) holds, if we replace Q k s with Pfe's, where Pfc is the fcth Legendre polynomial defined in 
Section |2~2"1 In other words, 



L, 



E4 n) ^ 



,/c=0 



00 

£ [( X „ - fc(fc + 1)) ajf > - c 2 (a[%A k + a^B k + a[%C k 

k=2 

+ [(Xn 1(1 + 1)) «i n) - C 2 + «i n) 5l) 



Pa 



( X „ - 0(0 + 1)) 4 n) - c 2 (4"> A, + 4">B ) 



P 
Po 



- c 2 (4 n) (t 2 P! - PxPi - C3P3) + 4 n) (t 2 p, - p P) - c 2 p 2 ; 

We combine ([75]). of Section to conclude that 

t 2 • P x (i) - P a • P x (t) - C 3 • P 3 (t) = 0, 
i 2 • P (i) - P ■ P (t) - C 2 • P 2 (t) = 0. 



(528) 



(529) 



We recall that {P k } form an orthogonal system in L 2 [—1,1], and combine this observation with (ptB"]) 
in Section O & in Section (|5M1) . (|52"5|i and (152T))) to conclude that, for every integer k > 2, 



( Xri - fc(fc + 1)) a™ - c 2 (ajg^ + 4" } P fc + atl 2 C k ) = 0, 



(530) 



and also 



( Xn - 1(1 + 1)) 4" } " C 2 (4"^1 + «i n) 5l) 

( X „ - 0(0 + 1)) 4 n) - c 2 (4"> A, + 4 n) P 



-0, 
= 0. 



(531) 
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We substitute (j53"0) , (j53T) into (jS"2"T) and use ([5T5)) to obtain 



L,, 



(t) =L n 



,fc=0 



(t) 



c 2 a^ ] (t 2 Qi(t) - B x Qx{t) - C 3 Q 3 (t)) 
c 2 a { n) (t 2 Q (t) - B Q Qo(t) - C 2 Q 2 (t)) 



(532) 



We combine (JHSJ), dHZl), ([TOO) - (fTUTl) of Section l2~2l to obtain 



t 2 Qo(t) - B Q (t) - C 2 Q 2 (t) 



and 



3/2 8 l-t 



t 2 Q!(t) - BxQxit) - C 3 Q 3 (t) 



(533) 



t . 1 + f 
- log 

2 6 1 - t 



5 V 4 



(5t 3 - 3t) log - 



l + t 



1 

To 



3*) J logy 



l+t , 3 o 4 

t 2 + - + t 2 

t 5 15 



- t 

a 



5 , 2 
-t 2 + - 
2 3 



(534) 



Finally, we substitute (j53"3) . (|5Mj) into ([5321 to obtain ([523)1 . 



In the following corollary, we establish a recurrence relation between the derivatives of $ n of 
arbitrary order (compare to Theorem 1151 in Section [2.1[) . 



Corollary 6. Suppose that the function $ n : (—1,1) 
Suppose also that — 1 < t < 1 is a real number. Then, 



(l-t 2 )-^(t)-4t-^(t) + ( Xn -c 2 t 2 



y R is defined via (|518[) o/ Theorem [FT} 
2) -<(t)-2c 2 t .* B (t) = 

(535) 



where is defined via 



Section \2.2\ (comvare to (]72|) o/ Theorem ] 15\ in Section \2.1\) . Also, 



for every integer k > 2, 

(l-t 2 ) $i fe + 2 ) (t) - 2 (fc + 1) t*<* +1 > (t) + ( x „ - k (k + 1) - c 2 t 2 ) $W (t) 

- cMj?- 1 ^*) - c 2 fc(fc - 1) $<*" 2) (t) = 0. (536) 

In other words, the higher order derivatives of $„ and ip n satisfy the same recurrence relation (|73[) 
(see Theorem \15\ in Section \2.1\) . 

Proof. To prove (|535[) . we differentiate both sides of (|525|) with respect to t. To prove ([536D . we 
observe that the second derivative of the right-hand side of (1525[) is identically zero, and combine 
this observation with Theorem [15] in Section 12.11 ■ 



The rest of this subsection is devoted to establishing the positivity of the quadrature weights 
Wi, . . . , W n , defined via (|434p in Definition [21 The principal result of this part is Theorem [731 (see 
also Remarks 112 1113) . 
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Theorem 69. Suppose that c > is a real number, and that n > is an odd integer. Suppose 
also that t\, t%, ■ ■ ■ , t n and W\, W2, • ■ • , W n are defined, respectively, via (|432l) , (I434[) in Definition^ 
Suppose furthermore that the integer jo is defined via the formula 



Jo = Z±±. (537) 



Then, for every integer j = 1, . . . , n, 



MM? Wi = W * + WT)Jo Mt)dt (538) 



Proof. Suppose that the differential operator L n is defined via (|524[) in Theorem [68J Suppose also 
that the function <& n : (— 1, 1) — » R is the solution of the homogeneous second-order ODE 

L n [<p] = (539) 
in the interval (— 1, 1) with the initial conditions 

$ " (0) = ^7o)' $ « (0) = - (540) 

Obviously, $ n is an even function. Moreover, 

$ n (t)-tf n {t)-& n (t).Mt) = Y=p (541) 

for all real — 1 < t < 1 (this is the classical Abel's formula; see e.g. Theorem 3.3.2 in [S]). Suppose 
that the function l>„ : (— 1, 1) -> M is defined via (|518[) in Theorem |57] We combine (|525|) of 
Theorem [SS] with (|598l) to conclude that <£„ satisfies the non- homogeneous ODE 



Ln 



(x) = , (542) 



for all real — 1 < x < 1. We observe that ip n ,Q n are two independent solutions of the ODE (|539l) . 
and combine this observation with (|542[) to conclude that, for all real — 1 < x < 1, 

*„(a:) = Ci ■ ^ n (a) + C 2 • $ n (a?) + 

"' A "'-'' (,)l ^ n (x) f X $ n (t)dt-<Z> n (x) [ X Mt)dt), (543) 



for some constants C\,C%. Out of the four summands on the right-hand side of (|543l) , the function 
C\ ■ ip n {x) is odd, while the other three functions are even. We combine this observation with (|518j) 



and (|543[) to conclude that 

Cx = 0. (544) 
On the other hand, we substitute x — into (|543[) to conclude that 

C2 = W (M5) 
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Suppose now that j is an integer between 1 and n. We recall that tj is a root of ipn due to (|432[) . 
and combine this observation with (|540p , (|543l) , (|544l) , (15451) to obtain 

= • (*»(0)<(0) - lcXn f (0) ^ tf n (t) di) . (546) 

We combine (|52"51) of Corollary [3] with (|537)) and ([516]) to obtain 

Wj ■ fl/ n (tj) = • fa, • (C(0)) 2 + icA»^(0) tf n (t) dt\ . (547) 

Finally, we combine ([BUI with (IMTf to obtain {535]). ■ 
Theorem 70. Suppose that c > is a positive real number, and that 

c > 30. (548) 
Suppose also that n > is an odd positive integer, and that 

2c 

n > — + 7. (549) 

Suppose also that t%, . . . , t n and W\, . . . , W n are defined, respectively, via (|432[) . (|434[) of Definition^ 
Suppose, in addition, that 

W {n+1)/2 <2-\\ n \-V2n~. (550) 

Then, 

W 1 + --- + W n < 4 ^'^") 7/4 ,| An |. (551) 

Proo/. We combine (|549[) . Theorems HI H71 in Section |2~T1 and (|462j) in the proof of Theorem IBT1 to 
conclude that 

Bjjiiji^ <552) 

We combine (|552l) with Theorem Q] of Section |2~T1 to conclude that, for any — 1 < x < 1, 

777^/ dt < 4|A„|^ / |V„(t)|dt<2.|Xn|.V2S. (553) 

WnW Jo JO 

We combine (I550p . (|553[) with Theorem |6"91 to conclude that, for every integer 1 < j < n, 

w< MM) 4-|AJ-V2n. (554) 

We combine (|549p . (|554p and Theorems HI [TTl in Section I2TT1 to conclude that, for every integer 
1 < j < n, 

Wj < _ X " 2 • 4 • |A„| • V2^. (555) 

Xn C ■ tj 

We combine (I549p with Theorems H] in Section |2~T1 to obtain the inequality 

n < Vx^- (556) 



Now (I55ip follows from the combination of (|555l) and (|556[) . 
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Theorem 71. Suppose that c > is a positive real number, and that 

c > 30. 

Suppose also that n > is a positive integer, and that 

2c - 
n > h 7. 

7T 



(557) 



(558) 



Suppose also that t\, . . . , t n and W%, . . . , W n are defined, respectively, via (|432[) . (|434[) of Definition® 
Then, 



H i ir„ ■ 2 - |A„ ■ ( 24 - log + 130 • ^e; 



(559) 



Proof. Suppose that the function I(t) : (—1, 1) —> K is defined via (|4f 71) in Theorem [551 Then, 



1 51 



i ^n(*i)'(*-*i) 



+ /(t)-Vn(*), 



(560) 



for all real —1 < t < 1. We integrate (I560p over the interval (— 1, 1) and use Theorem[T]in Section |2~T1 
Theorems [46l [56l [57] and Definition [2] to obtain 



(561) 



Wi + ■ ■ ■ + W n > 2 - |A„| • ( 24 • log — + 13^fe + 40c|A„| + 2^2. 

|A n | 



We combine (g2J|, (O, Theorem g] in Section O with ({555)1 to obtain 

40c|A„| < 40 V^hTc < 40V2^ ■ 

We combine ([557] . ([558")l . ([51)2"]) with Theorem H in Section O to obtain 

13^ + 40c|A„| + 2\/2< 130^. 

Now we substitute (|51)3"f into (|56T|) to obtain (|559"|) . 

Theorem 72. Suppose that c > is a positive real number, and that 

c > 30. 

Suppose also that the real number f3 is defined via the formula 

90 



/3 = 



log(30) ■ 

Suppose furthermore that n > is a positive integer, and that 

2c -log(c) , / 16ec 



The 



| A»l- f24.log -L + 130^T+ 4 ^ (X " ) 2 7/4 > | <2-e~ w . 



(562) 

(563) 
■ 

(564) 
(565) 

(566) 
(567) 
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Proof. We combine (15641) , (I565[) , (|566[) with Theorem @] in Section 12.11 to obtain the inequality 

4^( Xra ) 7 / 4 13Q.( X ») 5 / 4 + 4^( X „) 7 / 4 . 10-(x») 7/4 

loUyXn -f" „ <. „ <. 9 ■ 

Xn - C* Xn~ c Xn - C Z 

Also we combine (|564l) . (|565l) . (|566l) with Theorem [7] in Section |2~T1 to conclude that 

Xn > c 2 + P ■ log(c) • c. 
We combine (|5M1) . ([51)51 and to obtain 

(Xnf 4 < C 3 / 2 -(l+/3.l0g(c)/c) 3/4 



/3 • log(c) • c 



£ • log(c) ' 



We substitute ([5T0")) into ([51)81 to obtain 



130^ + 4 ^K /4 < 10 ^- X - 



p ■ log(c) 



We combine ([553]) . ([555]) . dSJgJ with Theorem [TT] to obtain 



Y 4 
An 



|A„| < 1195 - ^ - exp 



7T Xn ~ C 
4 ' 



'Xn 



(568) 
(569) 

(570) 
(571) 
(572) 



We combine ([5001) . ([gOTjl in the proof of Theorem [6j with (plj) . ([565]) . ([569"]) . (f57T]) . ((572)) to obtain 
( 4V2(Xn) ?/4 ^ 



|A„| • 130^fc + 



Xr, 



< 



11950 ■ cVgc ■ (1 + P ■ log(c)/c)" 
P ■ log(c) 

11950 -^VSc^ 5 

• exp 



• exp 



/8 • log(c) 



4 ^1 + p ■ log(c)/c 



< 



/? • log(c) 



7T_/3j lQg(c) 

8 



(573) 



We take the logarithm of the right-hand side of ([5731) and use (|564[) . (|565[) to obtain 



log 



log 



11950- c 3 ^- 4 



P ■ log(c) 

'ii950-\/8-4 5N 



• exp 



TT ■ P- l0g(c) 



/3-io g (c) y 

We combine (|573|) with (|574[) to conclude that 

|A„| • ( 130^a~ + 



7 it- P 
2 8~~ 



4V2(Xn) 7/r 

Xn - C 2 



8 



log(c) < -10. 



< e 



We combine (15M]) . ([5rJ5|) . ([575]) to conclude that 

|A„| < e" 16 . 

It follows from (J576j that 

24 • |An| • log -L < 24 • 16 • e- 16 < e" 10 . 
|A„| 

Now (I567P follows from the combination of (|575l) and (|577[) . 



(574) 

(575) 

(576) 
(577) 
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Theorem 73. Suppose that c > is a positive real number, and that 

c > 30. 

Suppose also that n > is a positive odd integer, and that 

2c 



2c f C \ 

n > — + 5 • log(c) ■ log (-J 



(578) 



(579) 



Suppose furthermore that W\ , . . . , W n are defined, via (|434j) of Definition Then, for all integer 
j = 1, . . . ,n, 



W j > 0. 



(580) 



Proof. Suppose first, by contradiction, that 

W(n+l)/2 <2-|A„|-\/2^. 
Then we combine ([575)1 . (15791) . <t55T|) with Theorems EDJ [71] to conclude that 

V2-(x™) 7/4 



(581) 



|A„| >Wl 



>2-|A„|- (24.1og^ + 130-^ 



in contradiction to Theorem 1721 Therefore, 

W (B+ i )/2 >2.|A n |-V^n. 



(582) 



(583) 



We combine (|583[) with Theorem [5§] and (|553[) in the proof of Theorem [70] to obtain, for every 

j = 1, . . . ,71, 



(Cfe)f ■(!-*?) 



icA r , 



(C(o)r 



>2-|A„ 



'2n- 



C(0) Jo 

cAr, 



>o, 



(584) 



where t\ , . . . , t n are defined via (|432l) in Definition^ Now (|580p follows directly from the combination 
of ([553)1 and ([55511 in the proof of Theorem [701 ■ 

Remark 12. The conclusion of Theorem \ 73\ holds for even integers n as well. The proof of this fact 
is similar to that of Theorem \ 73\ and is based on Theorems \ 71\ \72\ and the obvious modifications of 
TheoremWj\7% 

Remark 13. Extensive numerical experiments (see e.g. Table Wb\ and Figure [X^j) seem to indi- 
cate that the assumption (|579p is unnecessary. In other words, the weights W\ W n are always 
positive, even for small values of n. 

Remark 14. It follows from Theorem \69\ that, if 1 < j,k < n are integers, then 



«(*i)) 2 ■ (1 - £) • Wi = «(* fc )f ' (1 - 4) -W k + (|A„|) 



(585) 



78 



(see also Experiment 15 in Section \6.2.2\) . We observe that for c = the quadrature introduced 
in Definition [H is the well known Gaussian quadrature, whose nodes are the roots t\ , . . . , t n of the 
Legendre polynomial P n (see Section WJfy , and whose weights are defined via the formula 



(see e.g. \Tj, Section 25.4). Thus, (|585l) is not surprising. 

5 Numerical Algorithms 

In this section, we describe several numerical algorithms for the evaluation of the PSWFs, some 
related quantities, and the nodes and weights of the quadrature, defined in Definition[2]in Section l4~4l 
Throughout this section, the band limit c > is a real number, and the prolate index n > is a 
non-negative integer. 

5.1 Evaluation of Xn and ip n (x), ip' n {x) for — 1 < x < 1 

The use of the expansion of if> n into a Legendre series (see (1521 in Section |2"T2"|) for the evaluation of 
ip n in the interval (—1,1) goes back at least to the classical Bouwkamp algorithm (see [3]). More 
specifically, the coefficients , /3± , . . . of the Legendre expansion are precomputed first (see ((551) , 
(|84p in Section |2~2|) . These coefficients decay superalgebraically; in particular, relatively few terms 
of the infinite sum (|82l) are required to evaluate if) n to essentially machine precision (see Section [2~2l 
in particular Theorem [T5I and Remark [21 and also [55] for more details). 

Suppose that n > 0, and we are interested to evaluate the coefficients /3q , Pi , . . . of the 
Legendre expansion of ip m , for every integer < m < n. This can be achieved by solving two 
N x N symmetric tridiagonal eigenproblems, where N is of order n (see Theorem [T5I and Remark [5] 
in Section [2.21 and also [38] for more details about this algorithm). In addition, this algorithm 
evaluates Xo> • • • >Xn- Once this precomputation is done, for every integer < m < n and for every 
real — 1 < x < 1, we can evaluate ip m {x) in 0(n) operations, by computing the sum (1821) . 

Suppose, on the other hand, that we are interested in a single PSWF only (as opposed to all the 
first n PSWFs). Obviously, we can use the algorithm mentioned above; however, its cost is 0(n 2 ) 
operations (see Remark[2|). In the rest of this subsection, we describe an algorithm for the evaluation 

of /3g n \ (3[ n \ . . . and Xm whose cost is only 0{n) operations. 

This algorithm is also based on Theorem [T5J in Section 12.21 It consists of two principal steps. 
First, we compute a low- accuracy approximation x-n of Xn, by means of Sturm's bisection (see 
Section [2.7.51 ([M)l . ([Ml) and Remark [5] in Section [2^21 and also [5]). Second, we compute x« and 
ft n \ defined via ^ in Section O by means of inverse power method (see Section 12.7.41 and also 
[37) . [7]). The inverse power method requires an initial approximation to both the eigenvalue and 
the eigenvector; for this purpose we use, respectively, Xn and a random vector of unit length. 

Below is a more detailed description of these two steps. 

Step 1 (initial approximation Xn of Xn)- Suppose that the infinite symmetric tridiagonal 
matrices A even and A odd are defined, respectively, via ([50]) . (j^TT) in Section |2"T21 Suppose also that 
A^ n > is the N x N upper left square sub-matrix of A even , if n is even, or of A odd , if n is odd. 
Comment. N is an integer of order n (see Remark [2]). The choice 



2 



(586) 




N = 1.1 ■ c + n+ 1000 



(587) 



is sufficient for all practical purposes. 



79 



• use Theorems 0] [5] and [5] in Section 12.11 to choose real numbers xo < yo such that 

xq < Xn < Vo- (588) 

Comment. For a more detailed discussion of lower and upper bounds on \ni see, f° r example, 
|25j . |26j . See also Remark [TBI below. 

• use Sturm's bisection (see Section 12.7.51) with initial values £0,2/0 to compute \ n . On each 
iteration of Sturm's bisection, the Sturm sequence (see Theorem 1241) is computed based on the 
matrix A^ n ' (see above). 

Comment. We only require that \ n be a low-order approximation to \n in the following 
sense: Xn is closer to \ n than to any Xk with k =/= n. 

Remark 15. The use of Sturm's bisection as a tool to compute the eigenvalues of a symmetric 
tridiagonal matrix goes back at least to in the context of PSWFs, it seems to appear first in J13f . 



The cost analysis of Step 1 relies on the following observation. This observation is based on 
TheoremsEllHEliniEllEliu Section [2TTl as well as on extensive numerical experiments and asymptotic 
expansions (see, for example, [38], [30], [36], [25], [26]). 

Observation 1. Suppose that n > is an integer. 

If < n < 2c/ w, then (e.g. it seems < \o < c ) 

Xn+l ~ Xn = 0(c). (589) 

If n > 2c/ n, then 

Xn+i-Xn = 0{n). (590) 

Remark 16. If < n < 2c/tt, then we combine Theorems [^J E| in Section \2.1\ to obtain 

n ■ (n + 1) < Xn < c 2 . (591) 

We combine (|589p , (|591[) and Corollary in Section \2. 7.5\ to conclude that, in this case, the cost of 
Step 1 is 0(n ■ log(c)) operations. If, on the other hand, n > 2c/tt, then we combine Theorems^ 
Corollary [H in Section \2. 7. 5\ and (|590[) to conclude that, in this case, the cost of Step 1 is 0(n) 
operations. 

Step 2 (evaluation of Xn and /3^ n ^). Suppose that Xn is an approximation to Xn, computed in 
Step 1 (in the sense that x.n is closer to Xn than to any other eigenvalue Xk)- Suppose also that 
N is that of Remark [5] in Section |2~21 (see also Step 1 above, and, in particular, (|587p ). and that 
fin) e r n jg defined via in Section O 

• generate a unit length random vector /3 € R*. 

Comment. We use x„ and /3 as initial approximations to the eigenvalue Xn and the corre- 
sponding eigenvector, respectively, for the inverse power method (see Section [2.7.4[1 . 

• conduct inverse power method iterations until Xn is evaluated to machine precision. The 
corresponding unit eigenvector is denoted by (3^ . 

Comment. Each iterations costs 0(n) operations, and only O(l) iterations are required (see 
Section [2.7. 41) . In practice, the number of iterations is always less than 10. 
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conduct additional K iterations of inverse power method, until the convergence of the first 
coordinate of (3^ n \ 

Comment. Both analysis and numerical experiments (to be reported at a later date) suggest 
that 



K = 1 + ceil 




(592) 



where e is the machine precision (e.g. e w 1D-16 for double precision calculations), and ceil(a) 
is the minimal integer greater than a, for a real number a. For example, if |/3q I ~ 1D-99, and 
e « 1D-16, then K = 8. In practice, K does not to be known in advance; rather, we iterate 
until convergence. 

Remark 17. The cost of Step 2 is 0(n) operations. 

Remark 18. It is a well known fact (see e.g. \S7\ , J^j) that Xn is evaluated to essentially machine 
precision by the inverse power method. In other words, suppose that e is the machine accuracy 
(e.g. e ~ 1D-16 for double precision calculations); then, Xn is evaluated with relative accuracy 
e. In addition, (3^ approximates (3^ with relative accuracy e. However, this means that a single 
coordinate of is only guaranteed to be evaluated with absolute accuracy e. More specifically, 
for every integer k = 0, . . . , N , 



ft? 



in) 



< 



4 n) 



(593) 



We make the following observation from Remark [TSJ If we use (3^ to evaluate Legendre series 
(see ([521 in Section [2^21 and also (I595p . (|596[) below), the result will be obtained with high accuracy. 
On the other hand, the small coordinates of /3'"^ are only guaranteed to be computed with low 
accuracy. In particular, due to f|593[) . if, for example, |/3^,"' ) | < e/10 for some k, then, apriori, we do 
not expect to coincide with (3^ in any digit at all! 

The following conjecture states that the situation is much better than Remark [15] seems to 
suggest. This conjecture has been confirmed by both some preliminary analysis (see e.g. |27j . |28j ) 
and extensive numerical experiments. The matter is a subject of ongoing research; the results and 
proofs will be published at a later date. 

Conjecture 1. The coordinates of (3^ are evaluated with high relative accuracy. More specifi- 
cally, for every 1 < k < N , 



< e ■ log(VH), (594) 



where (3^ is the numerical approximation to (3^ n \ computed in Step 2, and e is the machine accuracy 
(e.g. e ~ 1D-16 for double precision calculations). 

In particular, Conjecture Q] implies that, no matter how small /3i™^ is, it coincides with (3^ in 
all but the last log 10 (\/c) decimal digits. 
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Evaluation of ip n (x), ip' n {x) for — 1 < x < 1, given \n and Pq U \ p[ n \ . . . Suppose Xn and 
the coefficients /3q , ... of the Legendre expansion of defined via ([85jl . in Section l2~2l have 
already been evaluated. Suppose also, that the integer AT is that of Steps 1,2 above (see, for example, 

H5S3)- 

For any real — 1 < x < 1, evaluate ^> n {x) via the formula 

^„ (x) = E (x) ■ a ( k n) = E ^ ( x ) ■ ■ VkTTj2. (595) 

k=0 k=0 

Also, we evaluate ip' n {x) via the formula 

2N 2N 

= E ^ • 4^ = E • ^ ■ + (596) 

fe=l fe=0 

Remark 19. The cost of the evaluation of \n o,nd /3q , /Sj , . . . wa Steps 1,2 is 0(n) operations 
(see Remarks 1161 117| above). Once this precomputation has been done, the cost of each subsequent 
evaluation ofip n (x), ip' n (x), for any real —1 < x < 1, is 0{n) operations, according to (|595[) . (|596|) 
and Remark^ in Section[ 



5.2 Evaluation of A n 

Suppose that the coefficients (3^ , f3^ , . . . of the Legendre expansion of ip n (see ([55)1 in Section 
as well as ip n (0), ip' n (0) have already been evaluated by the algorithm of Section I5TT1 If n is even, 
we compute A„ via the formula 

A » = -4^ / Mt) dt = ^ = ^4=-. (597) 



MO) J-l V-n(O) Vn(0) 

If n is odd, we compute A„ via the formula 

ic r 1 2 ica {n) l~2 icB [n) 

^ffi/j-^^^a'frffi (598) 

(see ([37]) in Section O and f75 ]) . ([50 ]) . ([55 ]) . ([51 ]) in Section EH]). 

Observation. According to (|597[) . (|598[) . the eigenvalue A„ is evaluated in O(l) operations as a 
by-product of Steps 1,2 of the algorithm of Section I5TT1 (the cost of these steps is 0(n) operations, 
due to Remarks [TBI [TT)) . Obviously, A„ and /3q , /8j are evaluated to the same relative accuracy. 
In particular, even though |A n | can be extremely small, A„ is evaluated with fairly high precision 
(see Conjecture [TJ in Section l5"7i"]) . 



5.3 Evaluation of the Quadrature Nodes 

Due to Definition [2] in Section l4~4l the n quadrature nodes t±, . . . , t n are precisely the roots of if> n in 
(—1,1). In this subsection, we describe a numerical algorithm for the evaluation of the quadrature 
nodes. Since tp n is symmetric about the origin (see Theorem[JJin Section I2TT]) . it suffices to evaluate 
the roots of i/j n in the interval (0, 1). 

To evaluate the quadrature nodes, we use the fast algorithm for the calculation of the roots of spe- 
cial functions, described in [UJ. This algorithm is based on Priifer transformation (see Section l2~6l) . 
Runge-Kutta method (see Section T2.7.3P and Taylor's method (see Section f2.7.2l) . It computes all 
the roots of ip n in (—1, 1) in only 0{n) operations. 
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A short outline of the principal steps of the algorithm is provided below. For a more detailed 

description of the algorithm and its properties, the reader is referred to 

The following observation is a direct consequence of Theorem [33J in Section 12.61 

Observation 1. Suppose that the function 9 : [ti,t n ] — > R is defined via (|148[) in Theorem |2"31 

in Section l2~6l Suppose also that the function s : [ir/2,ir ■ (n — 1/2)] — > [— t n ,t n ] is the inverse of 9. 

Then, s is well defined, monotonically increasing and continuously differentiable. Moreover, for all 

real tt/2 < rj < n ■ (n - 1/2), 

S ' {V) = /( s (^))- u ( s ( J7 )).sin(27 7 )' (599) 

where the functions f,v are defined, respectively, via (|144j) . (|145[) in Section |2"1H In addition, for 
every integer i = 1, . . . , n, 

i - A = *i. ( 6 °°) 



2 



and also 



-(t)=°- (601) 



Suppose now that i m j n is the minimal root of tp n in [0, 1). 



Step 1 (evaluation of t m in)- If"- is odd, then 

tmin = *(n+X)/2 = 0, (602) 

and this step of the algorithm is trivial. On the other hand, if n is even, we observe that 

tmin = t(n+2)/2 > 0. (603) 

We numerically solve the ODE (|599[) with the initial condition (|601l) in the interval [im/2, w ■ (n + l)/2], 
by using 20 steps of Runge-Kutta method (see Section l2.7.3[) . The rightmost value t m i n of the solution 
is a low-order approximation of i m ; n (see (|600l) . (|603I) V 

We compute t m ; n via Newton's method (see Section ^. 7.1[) . using t m i n as the initial approximation 
to t m in- On each Newton iteration, we evaluate ip n and tp' n by using the algorithm of Section [5. II 

Observation 2. The point t m ; n approximates t m [ n to roughly three-four decimal digits. Subse- 
quently, only several Newton iterations are required to obtain t m i n to essentially machine precision 
(see [11] for more details). Thus, the cost of Step 1 is 0(n) operations. 



Step 2 (evaluation of ip' n (t m i n )) . We evaluate V4(tmin) by using the algorithm of Section l5Jl 
Observation 3. The cost of Step 2 is 0(n) operations (see Remark IT91 in Section UTTj) . 
The remaining roots of i\) n in (i m i n , 1) are computed iteratively, as follows. Suppose that n/2 < 

j < n is an integer, and both tj and V4(ti) have already been evaluated. 

Step 3 (evaluation of and ip' n (tj+x), given tj and ijj' n (tj)). 

• use the recurrence relation (1731 (see Theorem[15]in Section l2~Tj) to evaluate ijjn (tj)> • ■ ■ > tyn (tj)- 

• use 20 steps of Runge-Kutta method (see Section 12.7. 3[) . to solve the ODE (|599|) with the 
initial condition 




(604) 
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in the interval {it ■ (j — 1/2), tt ■ (j + 1/2)], by using 20 steps of Runge-Kutta method (see 
Section [2.7.3p . The rightmost value tj+i of the solution is a low-order approximation of tj+i. 

• compute tj+i via Newton's method (see Section [2 .7. II) . using tj+\ as the initial approximation 
to tj+i. On each Newton iteration, we evaluate ip n and V'n by using Taylor's method (see 
Section f2. 7. 21) . The Taylor expansion of order 30 about tj is used, e.g. 

30 ,(k) ( , 

M t) = £ ■ (* - W k + 0((t- t,) fc+1 ) • (605) 

fe=0 

• evaluate V^fe'+i) by using Newton's method, i.e. by computing the sum 

29 ,(k+l), s 

E^r^'^-^- (606) 

k=0 

Observation 4. The point approximates tj+\ to roughly three-four decimal digits. Subse- 
quently, only several Newton iterations are required to obtain tj+i to essentially machine precision 
(see [TT] for more details). The cost of Step 3 is 0(1) operations. 

Step 4 (evaluation of tj and ip'n(tj) f° r a U j < n/2). Step 3 is repeated iteratively, for every 
integer n/2 < j < n. To evaluate tj and ipnftj) f° r — 1 < tj < 0, we use the symmetry of ip n about 
zero, established in Theorem[JJin Section HOI More specifically, for every 1 < j < n/2, we compute 

tj = t n+1 -j (607) 

and 

<(ti) = (-l) n+1 -i(WH). (608) 

Summary (evaluation of tj and ip' n (tj), for all j = 1, . . . , n). To summarize, to evaluate the 
roots of ifj n in (—1,1) as well as ip' n at these roots, we proceed as follows. 

• run Step 1, to evaluate t min (see (|602[) . (|603p ). Cost: 0(n). 

• run Step 2, to evaluate V'n^min)- Cost: 0(n). 

• for every integer n/2 < j < n, run Step 3. Cost: 0(n). 

• for every integer 1 < j < n/2, run Step 4. Cost: 0(n). 

Remark 20. We observe that the algorithm of this subsection not only computes the roots ti, . . . ,t n 
of ip n in (—1,1), but also evaluates ip' n at all these roots. The total cost of the algorithm is O(n) 
operations. 

5.4 Evaluation of the Quadrature Weights 

In this subsection, we describe an algorithm for the evaluation of the weights W\ , . . . , W n of the 
quadrature, defined in Definition [5] in Section 14.41 The results of this subsection are illustrated in 
Table l20l and in Figure [12] (see Experiment 15 in Section IB. 2. 2p . 

Obviously, one way to compute W\ , . . . , W n is to evaluate the integrals of ip\ , . . . , <p n numerically 
(see Definition [2]) . However, each tpj has n — 1 zeros in (—1, 1), and this approach is unlikely to 
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cost less that 0{n 2 ) operations (see also Section I3~T1) . In addition, each <pj has a singularity (albeit, 
removable) at tj , which might be a nuisance for numerical integration, especially if high precision is 
required. 

Below we describe two additional ways to evaluate the weights, based on the results of Sec- 
tion 14.4.41 One of them, based on Theorem [67] and Corollary [5] is straightforward and accurate; 
however, its cost is 0(n 2 ) operations. The other way, based on Theorem 1681 and Corollary [BJ in 
addition to having high accuracy and being easy to implement, is also computationally efficient: its 
cost is only 0(n) operations. 

We assume that the quadrature nodes t\, . . . , t n as well as ip' n (ti), . . . , tj)' n (t n ) have already been 
computed (by the algorithm of Section l5~3l whose cost is O(n) operations). 

Algorithm 1: evaluation of W\, . . . ,W n in 0(n 2 ) operations. Suppose that the coefficients 
dp , •••■> a 2N °f * ne Legendre expansion of ip n (see plj) in Section lz\2|) have already been evaluated, 
by the algorithm of Section 1531 here N is an integer of order n (see (|587[) in Section f5 . 1 1) . We 
compute Wj by evaluating the sum 

2N 

-^E^fe), (609) 

rnl J > k=0 

where Qo,Qi, - ■ ■ are the Legendre functions of the second kind, defined in Section l2~2l 

Observation 1. The sum (I609|) approximates the corresponding infinite sum to essentially 
machine precision, due to the superexponential decay of and the high precision to which Q/~ (tj ) 
are evaluated (see Sections 12. 2[ 15.11 and also [T2], [35], [T]). In combination with Theorem [57] 
and Corollary [S] this implies that (I609[) is an accurate formula for the evaluation of Wj (see also 
Experiment 15 in Section l"6.2.2[) . 

Observation 2. For every integer j, we evaluate Qo(tj), . . . , Q2N(tj) recursively, by using (|95l) . 
(IM1) in Section 12. 2[ in O(N) operations (see Remark [3] in Section 12. 2p . Since N — 0(n) (see 
Section [5TT]) . the overall cost of computing W\, . . . , W n via (I609P is 0(n 2 ) operations. 



Algorithm 2: evaluation of Wi, . . . ,W n in 0(n) operations. This algorithm consists of the 
following steps. 

Suppose that i m ; n is the minimal root of ij) n in [0, 1). In other words, 

imin = {^ 2 niS ° dd ' (610) 
[*(n+2)/2 n is even. 

Suppose also that the function $„ : (—1, 1) — > M is defined via f|518[) in Theorem l67l in Section 14.4.41 

Step 1 (evaluation of $„(t m in) and $^(t min )). Suppose that the coefficients a ,...,0.2$ of 
the Legendre expansion of ip n (see (|84p in Section l2"T2"j) have already been evaluated by the algorithm 
of Section [5TT1 We evaluate 3>n(^min) by computing the sum 

2N 

^a^OfeCW)- (611) 
fc=0 

Also, we evaluate $^(t m i n ) by computing the sum 

2N 

E4 n) QUw) (612) 



k=0 
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(see Algorithm 1 and Observations 1, 2 above, Theorem 1571 in Section [4.4.41 and Section [ 

Observation 3. We evaluate Q' (t min ), . . . ,Q' 2N (t m i n ) recursively (see Sections 12.21 15.11 and 

also [12], [38], [1]). Thus both (|611l) and (1612[) approximate $ ra (i m in) and ^' n (t min ), respectively, to 

essentially machine precision, and are computed in O(n) operations (see Observations 1, 2 above 

and Remark [3] in Section [2~2l . 

We evaluate $„ at all but the last four remaining roots of ip n in [0, 1) iteratively, as follows. 

Suppose that n/2 < j < n is an integer, and both $> n (tj) and & n (tj) have already been evaluated. 

Step 2 (evaluation of $ n (tj+x) and <b' n (t J+1 ), given & n (tj) an d 

• use the recurrence relation ([5351) . ([536D (see Corollary [B] in Section l4.4.4l) to evaluate $„ (t, ),..., $ 

• evaluate $ n (tj + i) by using Newton's method, i.e. by computing the sum 

60 ,(fc), - 

E^fcp-(^-*i) fe - ( 613 ) 

fc=0 

• evaluate $^(t J+ i) by using Newton's method, i.e. by computing the sum 

E^V^^ 1 "^ (614) 

fc=0 

Observation 4. For each j, the cost of the evaluation of (161 3[) . (I614[) is 0(1) operations (i.e. 
does not depend on n). Also, ([613[) . ([6141) approximate, § n {tj) and $^(£j), respectively, to essentially 
machine precision. For a detailed discussion of the accuracy and stability of this step, the reader is 
referred to [IT] . 

Step 3 (evaluation of $> n (tj) f° r n — 3 < j < n). For j = n — 3, n — 2, n — 1, n, we evaluate 
<& n {tj) by computing the sum 

E4 n) «»fe). (615) 

fc=0 

(similar to ([61ip in Step 1). 

Remark 21. We compute $ n at the last four nodes via ([615[) rather than (I613p . since i/ie accuracy 
of the latter deteriorates when (1 — n) becomes too small (see ([536[) in Corollary^. Since this 
approach works in practice, is cheap in terms of the number of operations and eliminates the above 
concern, there was no need in a detailed analysis of the issue (see also 111] for more details). 



Step 4 (evaluation of i& n (tj) for 1 < j < n/2). Suppose that 1 < j < n/2. We evaluate § n {tj) 
via the formula 

= (-1)" +1 • $„(tft +1 -i) (616) 

($ is symmetric with respect to zero due to the combination of Theorem [57J in Section 14.4.41 and 
P>]1 in Section [2~2l. 
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Step 5 (evaluation of Wi, . . . , W n ). By performing Steps 1-4 of Algorithm 2, we evaluate <t n 
at the roots t\,...,t n of -0n in (— 1, 1). Now, for every j — 1, . . . , n, we evaluate via (|523|) of 
Corollary [5] in Section |4"X41 

Remark 22. TTie overall cost of Steps 1-5 of Algorithm 2 is 0{n) operations. 
5.5 Evaluation of ip n and its roots outside (—1, 1) 

The PSFWs provide a natural way to represent bandlimited functions over the interval (—1, 1) (see 
Theorem [T] in Section l2.ip . Therefore, even though each ip n is defined (and holomorphic) in the 
whole complex plane, in applications (construction of PSWFs, quadratures, interpolation etc.) one 
is mostly interested in the properties of ip n {t) for real t inside (—1, 1) (see, for example, Section |2~T1 

EH], ES], EB], [23, 

On the other hand, the properties of the quadrature rules studied in this paper (see Definition [5] 
in Section [4.4[1 depend, perhaps surprisingly, on the behavior of ip n outside the interval (—1, 1) (see 
Sections 14.2.21 14.31 14. 4p . Thus, while one is rarely interested in the evaluation of ip n and related 
quantities outside (—1,1) per se, we do need such tools to illustrate our analysis (see Section [5] 
below) . 

The rest of this section is devoted to the description of numerical algorithms for the evaluation 
of ipn{x) and ^' n {x) for x > 1, as well as the location of the roots of -0n in (1, oo). These algorithms 
were developed as auxiliary tools, and are not meant to be used in practical applications. 

Throughout this subsection, we assume that c > is a positive real number, and n is a non- 
negative integer. 

5.5.1 Evaluation of tl> n (x) for x > 1 

To evaluate ip n {x) for x > 1, we use the integral equation (I37[) in Section 12.11 (as opposed to using 
the Legendre series ([52")) of Section [2^21 to evaluate ip n { x ) for —1 < x < 1). Namely, we evaluate 
ip n {x) via proceed as follows: 

• Compute Xn and the coefficients cXq , ol^, . . . of the Legendre expansion of tp n (see Sec- 
tion EH}. 

• Compute A„ (see Section IQ)) . 

• Compute ^ n {x) via evaluating the integral 

J- f Mt) ■ e 2CXt dt (617) 

*n J -I 

numerically, by using m = O(n) Gaussian quadrature nodes in the interval (—1, 1). 

We observe that the integrand in (|617[) is oscillatory: ip n has n zeros in (—1, 1), and e lcxt is periodic 
with period (2tt)/(cx). Moreover, ipn(x) itself is oscillatory with frequency of order n (unless x is 
between 1 and y'x^/c, see Theorems |2"51 1521 in Section l4.1.1j) . 

Thus, we used a fairly large number of Gaussian nodes to evaluate (|617p . For example, for 
c = 100 and n < 100 we used the Gaussian quadrature of order 500; for c = 1000 and n < 750 we 
used the Gaussian quadrature of order 3000. 

Remark 23. For each of the m Gaussian nodes t^, we compute ip n (Tk) via evaluating the sum 

2N 

J2 P M-4 n \ (618) 

j=0 
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where N is of order n (see Section [57J\) . Thus, the resulting algorithm for the evaluation of ip n (x) is 
fairly expensive: its cost is 0(N ■ n) = 0(n 2 ) operations, as opposed to 0(n) operations to evaluate 
ip n {x) for — 1 < x < 1 (see Remark[T$in Section [d^7]) . 



5.5.2 Evaluation of ip' n (x) for x > 1 

We differentiate the identity (|37p in Section HOI to obtain, for all complex x, 

CW = f / t ■ Mt) ■ e icxt dt. (619) 

We use (|619p to evaluate ip' n (x) for x > 1 in the same manner we use (|617|) to evaluate ip n (x) (see 
Section l5.5.ip . The resulting algorithm has the same cost as the one of Section l5.5.1l fsee Remark l23l) . 



5.5.3 Evaluation of the roots of ip n in (l,oo) 

Suppose that Xn > c 2 . Suppose also that k > 1 is an integer. According to Theorem |2"91 of 
Section mm 



= xq < xi < X2 < ■ ■ ■ < Xk- (620) 

V 

where x\, . . . ,Xk are the k minimal roots of ip n in (1, oo). We define the function 9 : [xq, Xk] R 
via (|22ip in Theorem [BUI of Section l4.1.1l Then, 9 is monotonically increasing; moreover, 

0(*o) = -|, ^i) = |, 6(x k )=ir-(k-±y (621) 

Also, 9 satisfies the nonlinear first order ODE (|222p (see Theorem I3TJ1) . Furthermore, for every 
integer j = 0, 1, . . . , k — 1, 

7T 

x, J+ i — Xj w — (622) 

(see Theorems [3ll [32] in Section |4~T1 for a more precise statement). 

Suppose now that j is an integer between and k— 1, and ceo, . . . ,Xj have already been evaluated 
(note that to evaluate the special point xo we only need to evaluate Xm see Section lOT) . We evaluate 
Xj+i as follows. 

• Define h via the formula 

h =ik- (623) 

• Use Runge-Kutta method (see Section [2.7. 31) to evaluate 9(xj + i ■ h) numerically (by solving 
the ODE ([2"2"2"|) with the initial condition (|B2"nO . for i = 1, 2, 3, ... . 

Comment. Due to (|622p . ft, defined via (1623)) is a reasonable step size of the Runge-Kutta 
ODE solver. 



Stop when 
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• Define ij+i via the formula 

x j+1 = x k + + ^ V h, (625) 
where i is as in (|624l) . This is the initial approximation of Xj+%. 

Comment. Due to (|622D . (|623D . we expect Xj+i to approximate Xj+i roughly to three-four 
decimal digits. 

• Use Newton's method (see Section 12. 7.1[) with the initial point Xj+i to evaluate Xj+%. 
Comment. For each Newton iteration, we evaluate ip n (x), ip' n (x) by using the algorithms of 
Sections 15.5.11 15.5.21 respectively. 

Remark 24. We observe that the algorithm of Section 1 5. 5. S\ is similar to that of Section 1 5.5*1 
However, rather than solving the ODE for the inverse of 9 (see (|599[) in Section ] 5. 3]) . here we solve 
the ODE for 9. Also, rather than evaluating ip n (x) and ipn( x ) by Taylor's method (see f|605[) . (|606[) 
in Section \5. 3\) . here we evaluate ip n (x) and ip' n (x) by using the algorithms of Section \5.5. 11 \5.5.2[ 
respectively. 

6 Numerical Results 

This section has two principal purposes. First, we illustrate the analysis of Section 0] by means 
of numerical examples. Second, we demonstrate the performance of the algorithms presented in 
Section [5j All the calculations were implemented in FORTRAN (the Lahey 95 LINUX version). 

In all the experiments, the principal numerical algorithms of the paper, described in Sections l5.1| - 
15.21 were run in double precision. On the other hand, the auxiliary algorithms of Section [53] (whose 
sole purpose is to illustrate the analysis) were run in extended precision. 

6.1 Properties of PSWFs 

In this subsection, we illustrate the analytical results from Section 14.11 Section 14.21 and Section 14.31 
6.1.1 Illustration of Results from Section 14.11 

Experiment 1. In this experiment, we illustrate Theorem [TBI in Section |2~TI and Theorem I2U1 in 
Section 14.1.11 We proceed as follows. We choose, more or less arbitrarily, the band limit c > and 
the prolate index n > 0, and evaluate ip n (x) at 1000 equispaced points in the interval (—1.5,1.5). 
To evaluate ipn{x) ioi —1 < x < 1, we use the algorithm of Section [531 (in double precision). To 
evaluate tp n (x) for \x\ > 1, we use the algorithm of Section [5.5.11 fin extended precision). 

We display the results of the experiment in Figures [TJ [5J corresponding to the choice c = 20, 
n — 9 and c = 20, n — 14, respectively. Each of these figures contains a plot of the corresponding 

We observe that the relations ([T4"|) and (f75|) hold for the functions in Figures [TJ [21 respectively. 
The inequality ([6"9")l of Theorem Q2] in Section 12.11 holds in both cases, that is, the absolute value 
of local extrema of ip n (t) increases as t grows from to 1. On the other hand, ([TO")) holds only for 
Figure[2j This is due to the fact that X9 < c 2 an d Xi4 > c 2 ( see also Theorem[4]in Section |2~T1) . Also, 
we observe that the magnitude of the oscillations outside (—1,1) is roughly inversely proportional 
to |A„|. 
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-1.5 -1 -0.5 0.5 1 1.5 



Figure 1: The function ip n (t) for c = 20 and n = 9. Since \n ~ 325.42 < c 2 , the behavior is as 
asserted in (|74[) of Theorem 1 1 61 The points ^Jx7i/c ~ 0.90197 and 1 are marked with asterisks. The 
eigenvalue |A n | 0.55978 is relatively large, and the oscillations of tp n outside (—1,1) have small 
magnitude. Compare to Figured Corresponds to Experiment 1. 



c 


n 




7T / x\ — 1 

2c V *i-(x»/c 2 ) 




10 
10 
10 


15 
19 

24 


0.46561E+00 
0.51090E+00 
0.55570E+00 


0.22542E+00 
0.24279E+00 
0.26055E+00 


0.20655E+01 
0.21043E+01 
0.21328E+01 


100 
100 
100 


76 
84 
92 


0.49260E-01 
0.57274E-01 
0.63570E-01 


0.23935E-01 
0.27070E-01 
0.29602E-01 


0.20581E+01 
0.21158E+01 
0.21475E+01 


1000 
1000 
1000 


652 
664 
676 


0.52819E-02 
0.56889E-02 
0.63367E-02 


0.23016E-02 
0.27295E-02 
0.30338E-02 


0.22949E+01 
0.20843E+01 
0.20887E+01 



Table 1: The relation between the left-hand side and the right-hand side of the inequality (|228|) 
of Theorem \31\ For each value of the band limit c, the three values of n are chosen such that 
|A n | ~ 10~ 5 , 10~ 9 , 10 -13 , respectively. Corresponds to Experiment 2. 
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v n (t) 




Figure 2: The function ipn{t) for c = 20 and n — 14. Since \n ~ 437.36 > c 2 , the behavior is 
as asserted in (|75[) of Theorem \16] The points 1 and yjxn/c ~ 1.0457 are marked with asterisks. 
Observe that \X n \ « 0.12564, and the oscillations of^ n outside (—1, 1) have relatively large magnitude 
(of order \X„\ ). Compare to Figure[J\ Corresponds to Experiment 1. 
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Experiment 2. In the following numerical experiment, we illustrate Theorem I3T1 in Section l4.1.1l 
We proceed as follows. For each of the three values of band limit c (namely, c = 10, 100, 1000), we 
pick three values of the prolate index n. The values of n are chosen to satisfy n > 2c/ n (which 
implies that \n > c 2 , due to Theorem H] in Section |2~TT) . Then, we evaluate the eigenvalue \n of the 
ODE (l48l) of Section |2~T1 by using the algorithm of Section IS~T1 Also, we evaluate the minimal root 
X\ of ip n in (l,oo) (see Theorem |2"91 in Section |4.1.1|) . by using the algorithm of Section f5. 5.31 

The results of this experiment are displayed in Table [T] This table has the following structure. 
The first two columns contain the band limit c and PSWF index n. The third column contains 
the difference between the x\ and the special point ^fx~nl c ( see Theorem l2"9l in Section B.l.ll) . This 
difference is the left-hand side of the inequality (|228D of Theorem [2D On the other hand, the fourth 
column contains the right-hand side of (|228|) (a lower bound on this difference). The last column 
contains the ratio of the value in the third column to the value in the fourth column. 

We observe that the value in the fourth column is smaller than the value in the third column 
roughly by a factor of 2, for all the choices of c, n. In other words, the lower bound on x\ — ^fxTi/c, 
provided by Theorem 1311 is rather inaccurate, but is of correct order. 
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Table 2: Illustration of Theorem \32\ with c = 100 and n = 90. |A n | w 10~ 10 . Corresponds to 
Experiment 3. 



Experiment 3. In the following numerical experiment, we illustrate Theorem I3"2l in Section Rl.l.ll 
We proceed as follows. We choose the band limit c and the prolate index n. For each such choice, we 
compute the first 20 roots x\, . . . , 220 of tp n in (1, 00), using the algorithm of Section l5.5.3l Also, for 
each k = 1, . . . , 19, we compute the upper and lower bound on Xk+i — Xk, established in Theorcm l32l 
The results of the experiment are displayed in Tables [2 [3l that correspond to c = 100 and 
n = 90, 110, respectively. These tables have the following structure. The first column contains the 
index k of the root Xk of ip n in (l,oo). The second column contains the difference between two 
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Table 3: Illustration of Theorem \ 3S\ with c — 100 and n = 110. |A„| ~ 10 25 . Corresponds to 
Experiment 3. 

consecutive roots Xk+i and Xk of ip n in (1, oo). The third and fourth columns contain, respectively, 
the lower and upper bound on Xk+i — Xk, as in (|233[) of Theorem 1321 The last two columns contain 
the relative errors of these bounds. 

We observe that the upper bound is more accurate in terms of relative error. Moreover, the 
relative accuracy of both bounds improves monotonically as k grows. On the other hand, for a fixed 
k, the accuracy in Table [H is slightly higher than that in Table [31 which suggests that the bounds 
worsen as n grows. We also observe that the difference Xk+i — Xk between two consecutive roots 
decreases monotonically to ir/c, as k grows (see (1235[) in Theorem 1321 and Remark [7]). 

6.1.2 Illustration of Results from Section 14.21 

Experiment 4. In this experiment, we illustrate Theorem [Ml m Section 14.2.11 We proceed as 
follows. We choose the band limit c = 10 and the prolate index n = 8. Then, we compute Xn by 
using the algorithm of Section 15.11 Also, we evaluate ip n and ip' n at 500 equispaced points in the 
interval 




(626) 



For each such point x, we compute Q(x) and Q(x), where the functions Q, Q are defined, respectively, 
via (gllD , ([M51 in Theorem IM1 

In Figures [3l [H we plot, respectively, Q and Q over the interval f|626|) . We observe that, as 
expected, Q is monotonically decreasing and Q is monotonically increasing. On the other hand, we 
observe that the second derivative of each of Q, Q does not have a constant sign in this interval. 
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Tabic 4: Illustration of Theorem \35\. with c — 100, n = 80, |A n | = 0.58925E-07. See Experiment 5. 
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Table 5: Illustration of Theorem \35l c = 200, n = 160, |A„| = 0.17136E-13. See Experiment 5. 
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Experiment 5. In the following experiment, we illustrate Theoreml551in Section l4.2.1l We proceed 
as follows. We choose, more or less arbitrarily, the band limit c and the prolate index n. For each 
choice of c, n, we evaluate Xn by using the algorithm of Section fSTTT Then, we evaluate the first 20 
roots x\, . . . , £20 °f m (1, °°), by using the algorithm of Section 15.5.31 (in extended precision). 
For each such root Xk, we evaluate ip' n (xk) by using the algorithm of Section 15.5.21 (in extended 
precision). 

We display the results of the experiment in Tables HI [SJ corresponding to c = 100, n — 80 and 
c = 200, n = 160, respectively. These tables have the following structure. The first column contains 
the index k of the root Xk of ip n in (1, 00). The second column contains the absolute value of the 
ratio of ip' n ( x k+i) to ip' n (xk)- The third and fourth columns contain the lower and upper bound 
on that ratio, respectively, established in (|254l) of Theorem [35] The last two columns contain the 
relative errors of these bounds. 

We observe that the ratio in the second column is always less than one. Moreover, it first 
decreases up to a certain k and then increases as k grows. Both bounds have roughly the same 
relative accuracy and become sharper as k grows. Even for k = 1 the errors are about 20%, while 
already at k = 7 they drop to about 5%. We also observe (not shown in the tables) that the 
magnitude of iV^^fc)! is about 10 8 for Tableland about 10 15 for Table[5](see also Experiment 6 
below) . 
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Table 6: Illustration of Theorem \41\ with c — 100, n = 100. A„ = 0.94419E-18. See Experiment 6. 



Experiment 6 In this experiment, we illustrate Theorems |4"T1 14"21 in Section |4.2. 21 We proceed as 
follows. We choose, more or less arbitrarily, the band limit c and the prolate index n. For each such 
choice, we evaluate \n and A„ , by using the algorithms of Sections 15.11 15.21 respectively (in double 
precision). Then, we compute the first 20 roots x\, . . . , X20 of ip n in (1, 00), by using the algorithm 
of Section 15.5.31 (in extended precision). For each such root Xk, we evaluate tp' n (xk) by using the 
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12750E- 
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23 
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23 
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23 





20502E- 


04 
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-01 


18 





13551E- 


23 





19164E- 


23 
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13745E- 
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19438E- 
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16479E- 
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10014E4-01 



Table 7: Illustration of Theorem\4j\with c = 1000, n = 700. X n = 0.12446E-21. See Experiment 6. 

algorithm of Section [5.5.21 (in extended precision). 

We display the results of the experiment in Tables H3 [3 corresponding to c = 100, n — 100 and 
c = 1000, n = 700, respectively. These tables have the following structure. The first column contains 
the index k of the root x k of ip n in (l,oo). The second column contains the reciprocal of |^(xfc)|. 
The third column contains the quantity 

ww- 1 ) 3 ' 4 (627) 



(4 - (Xn/C 3 ))" 1 |<Ml)|%/2 

(see (|306l) in Theorem l41"j) . The fourth column contains e k , defined via the formula 



2-K(l)l(**-(Wc 2 )) 



1/4 



\<(*k)\ = — '\\ % V " 4 " • (1 + e k ) (628) 
I An | (4-1) 37 

(we observe that e k in fl628[) is obtained via multiplying (|627[) by |^4(xfc)|/V^ and subtracting 
1 from the result). The last column contains b c {x k , y/Xn/c), defined via (|266[) of Definition [T] in 
Section 

According to Theorem 1411 the product of the values in the third and fifth columns is an upper 
bound on iV'nC^fc)! -1 (the second column). However, (|627[) alone (the third column) already over- 
estimates |V , n( JC *)l _1 by roughly v2- We also observe (see the fourth column) that the parameter 

defined via (1628[) . is fairly small, and decreases as k grows. According to Theorems l49l IST1 in 
Section 14.3.21 we expect e k to tend to zero as k grows to oo, since 

2-ivqi)i(4-(xn/c 2 )) 1/4 2.|^(i)i 

ttt ~ -p: r-, k y oo. (629) 

\\ n \(xl-lf 4 \K-x k \ 
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On the other hand, the fact that e k ~ 10 -4 already for k = 7 is somewhat surprising. In other 
words, the left hand side of (|629p is a fairly tight estimate of \ip' n (xk)\, even for small fc. 

We also observe that b c (xk, y/Xn/c) (see the last column) is very close to 1 even for fc = 1, and 
becomes even closer to 1 as fc increases. In other words, the upper bound e 1 / 4 ps 1.284 on this 
quantity (see Theorem in Section |4.2.2[) is somewhat overcautious. 



6.1.3 Illustration of Results from Section 



L,j=k 



j=k (i- Xj )^' n ( Xj ) 



,1/4 



y 



(z+iy dz 

( z 2- iXn/c 2 } f/2 



ratio 



1 

3 

5 

7 

9 

11 

13 

15 

17 

19 

21 

23 

25 

27 

29 

31 

33 

35 

37 

39 



0.29442E-19 
0.99172E-20 
0.57139E-20 
0.39054E-20 
0.29098E-20 
0.22851E-20 
0.18596E-20 
0.15530E-20 
0.13226E-20 
0.11441E-20 
0.10021E-20 
0.88694E-21 
0.79193E-21 
0.71242E-21 
0.64507E-21 
0.58742E-21 
0.53762E-21 
0.49425E-21 
0.45620E-21 
0.42261E-21 



0.31341E-17 
0.85727E-18 
0.46271E-18 
0.30749E-18 
0.22656E-18 
0.17760E-18 
0.14509E-18 
0.12209E-18 
0.10503E-18 
0.91920E-19 
0.81564E-19 
0.73193E-19 
0.66300E-19 
0.60534E-19 
0.55645E-19 
0.51451E-19 
0.47818E-19 
0.44643E-19 
0.41846E-19 
0.39365E-19 



0.10645E 
0.86442E 
0.80980E 
0.78735E 
0.77861E 
0.77720E 
0.78022E 
0.78614E+02 
0.79407E+02 
0.80345E 
0.81393E 
0.82524E 
0.83720E 
0.84969E 
0.86261E 
0.87588E 
0.88944E 
0.90325E 
0.91727E 



03 
02 
02 
02 
02 
02 
02 



-02 
-02 
-02 
-02 
-02 
-02 
-02 
-02 
-02 
-02 



0.93147E+02 



Table 8: Illustration of Theorem \44\ with c — n — 100. A n = 0.94419E-18. See Experiment 7. 



Experiment 7. In this experiment, we illustrate Theorem in Section 14.3.11 We proceed as 
follows. We choose the band limit and the prolate index to be, respectively, c = 100 and n = 100. 
We evaluate Xn and A„, by using the algorithms of Sections l5. 1115.21 respectively (in double precision). 
Then, we compute the first 40 roots x\ , . . . , X40 of ip n in (1, 00), by using the algorithm of Section r5.5.3l 
(in extended precision). For each such root Xk, we evaluate ip' n (xk) by using the algorithm of 
Section 15.5.21 (in extended precision) . 



For each fc = 1, 3, 5, ... , 39, we evaluate 

1 1 



max 

-KKl 



(t - x k ) ■ tp' n (x k ) (t - x k+1 ) ■ tp' n (x k+1 ) 



(630) 



(it turns out that the maximum is attained at t = 1.) Then, we evaluate the upper bound on (|630l) . 
provided by Theorem HH 

We display the results of the experiment in Table [5] The first column contains the index fc of the 
root Xk of if) n in (1, 00). The second column contains the quantity (I630p . The third column contains 
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the upper bound on (|630[) . provided by Q314I) in Theorem The last column contains the ratio of 
the third column to the second column. 

We observe that the quantity (|630[) (in the second column) decreases with k. On the other 
hand, the bound on (|630[) (in the third column) gets tighter as k increases from 1 to 11, and then 
deteriorates, as k increases further on, roughly linearly in k. The latter observation is not surprising, 
since 



|A„| 



(z + l) 2 dz 



7T ' |A„| 



I* (z 2 - (Xn/C 2 )) 3/2 C ' X* ' 

due to Theorem I3"2l in Section 14.1.11 while, for sufficiently large k 



k 



(t - X k ) ■ 1p' n (Xk) (t- X k +l) ■ 1p' n {x k+1 ) 



< 



20 ■ 7T • lAt, 



(631) 



(632) 



due to Theorem [52] in Section 14.3.21 In other words, the upper bound on (|630[) . provided by Theo- 
rem !44[ is of a wrong order (C^a;^ 1 ) instead of 0(x~^ 2 )). In particular, it can be used only to bound 
the head of the convergent series 



E 

k=l 



(t - X k ) ■ ip' n (xk) 



(633) 



Of course, this is precisely how Theorem |H] is used (see the proof of Theorem [55] in Section 14.3.11 
and the proof of Theorem [S3] in Section |4.3.3[) . 



8 
6 
4 
2 


-_ 2 i 



x 10 



■The whole series 



— First 50 terms 




-0.5 







0.5 



Figure 5: Illustration of Theorems El with c = 100, n = 80. |A„| = 0.58925E-07. 
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0.99408E-08 
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0.55502E- 


h03 


100 
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0.28195E-08 


0.19431E-07 
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0.58207E- 


-03 


100 


90 


0.63405E-13 


0.45487E-12 


0.71741E- 


fOl 


0.84186E- 


-03 


100 


91 


0.14648E-13 


0.12985E-12 


0.88645E- 


fOl 


0.87239E- 


-03 


200 


146 


0.57204E-08 


0.32856E-07 


0.57436E- 


fOl 


0.62129E- 


-03 


200 


147 


0.19902E-08 


0.12477E-07 


0.62691E- 


fOl 


0.64480E+03 


200 


158 


0.21537E-13 


0.15123E-12 


0.70219E- 


fOl 


0.91959E- 


-03 


200 


159 


0.64626E-14 


0.51123E-13 


0.79107E- 


fOl 


0.94591E- 


-03 


400 


274 


0.15108E-07 


0.80630E-07 


0.53369E- 


fOl 


0.67438E- 


-03 


400 


275 


0.61774E-08 


0.34713E-07 


0.56193E- 


fOl 


0.69478E- 


-03 


400 


288 


0.47053E-13 


0.31193E-12 


0.66293E 


fOl 


0.97598E- 


-03 


400 


289 


0.17000E-13 


0.12189E-12 


0.71703E- 


fOl 


0.99872E- 


-03 


800 


530 


0.18269E-07 


0.91984E-07 


0.50351E- 


fOl 


0.77801E+03 


800 


531 


0.83405E-08 


0.43433E-07 


0.52075E- 


fOl 


0.79612E- 


-03 


800 


546 


0.46822E-13 


0.29701E-12 


0.63434E+01 


0.10833E- 


-04 


800 


547 


0.19631E-13 


0.12945E-12 


0.65942E^ 


fOl 


0.11033E- 


-04 



Table 9: Illustration of Theorem 1 58[ See Experiment 8. 
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Experiment 8. In this experiment, we illustrate Theorem [55] in Section 14.3.31 We proceed as 
follows. We choose, more or less arbitrarily, the band limit c and the prolate index n. Then, we 
evaluate \n and A n , by using the algorithms of Section I5~T1 15.21 respectively (in double precision). 
Next, we find the roots t\, . . . , t n oiip n in the interval (—1, 1), by using the algorithm of Section I5~3l 
(in double precision). For each root ti, we compute ip' n (ti). 

Suppose now that the function / : [—1, 1] — > R is defined via (|417[) in Theorem [56] We evaluate 
I at 3 • (n + 1) points Z\, . . . , ^3( n +i) hi the interval [—1, 1]. The points are chosen in such a way 
that, if tk < Zj < tk+i for some j, k, then 

\ < Zj ~ tk < 3. (634) 
3 tfc+i — Zj 

In other words, no point Zj is "too close" to any root of ip n in (—1, 1). For each j = 1, . . . , 3 • (n+ 1), 
we evaluate I(zj) in extended precision. 

Remark 25. For each —1 < t < 1, we expect I(t) to be of order \X n \, due to Theorems \5b\ \5R On 
the other hand, suppose that —1 < t < 1, and tk is the closest root of ip n to t. Then, 

' ' 0(1). (635) 



Therefore, in the evaluation ofip n (t), we expect to lose roughly 

l0Sl ° (k(^U)-aJ (636) 

decimal digits. In other words, this calculation is rather inaccurate. However, since we need it only 
to illustrate the analysis, we were satisfied when we got at least two decimal digits, and did not make 
any attempts to enhance the accuracy. 

On the other hand, we compute the first 50 roots x\, . . . , X50 of ip n in (1, 00), and, for each such 
root Xj, we evaluate tp' n (xj). These calculations are based on the algorithms of Sections 15.5.21 15~5.3I 
Then, for each Zj, we evaluate the sum 

50 / 1 1 \ 

ho(zj) 2J \^' n ( Xk ) ■ (zj - x k ) + i>' n (-Xk) ■ (zj + x k ) ) ' ^ 37 ' ) 

We display the results of the experiment in Figures [S] El for c = 100, n = 90 and c = 100, n = 91, 
respectively. On each of these figures, we plot the function /, defined via (|417p in Theorem [56] (blue 
solid line) and the function /50, defined via (I637|) (red dashed line). 

We observe that, in both figures, the maximum of both / and /50 is attained at the end points 
of the interval. Also, we observe that the values of I and /50 are of order |A n |, as expected; also, 
the functions appear, at least by eye, to be well approximated by polynomials of order up to 3. In 
other words, the reciprocal of ip n seems to be approximated up to an error of order |A„| by a rational 
function with n poles, as asserted in Theorems 1561 l58l 

We display additional results of this experiment in Tabled This table has the following structure. 
The first and second column contain the band limit c and the prolate index n, respectively. The 
third column contains the maximum of the absolute value of the function / in the interval [—1, 1], 
i.e. 

||J|| 00 =max{|7(t)| : -!<*<!}, (638) 
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where / is defined via (I417P in Theorem [55] The fourth column contains | A„ | . The fifth column 
contains the ratio |A n |/||/||oo. The last column contains I maX) defined via (|419|l in Theorem 1561 

We make the following observations from Table GO First, |A n | alone is already an upper bound 
on ||/||oo- Moreover, for a fixed band limit c, the ratio |A n |/||/||oo increases as n grows. For all the 
values of c, n in Table [9j this ratio varies between 5 and 9. On the other hand, I max varies between 
500 and 1000. Moreover, I max increases with n, for each fixed band limit c. In other words, the 
upper bound |A n | • / max on ||/||oo, established in Theorem 1561 deteriorates as n increases. Moreover, 
the factor / max in (|418[) of Theorem [55] appears to be unnecessary. The main source of inaccuracy 
is Theorem in Section 14.3.11 which provides a relatively poor upper bound on the expressions of 
the form (|637[) (see Figures [5j [6] and Experiment 7 above) . 

Nevertheless, due to the fast decay of |A„| with n, the estimates of Theorem IS"6l albeit somewhat 
loose, are sufficient for the purposes of this paper (see the analysis of the quadrature error in 
Section 03] and also Experiment 14 in Section 16.2.11 below). 



6.1.4 Illustration of Results from Section 14.41 

Experiment 9. In this numerical experiment, we illustrate Theorem [SHI in Section 14.4.11 We 
proceed as follows. We choose, more or less arbitrarily, the band limit c, the prolate index n and 
the root index 1 < j < n. Then, we evaluate A„ and the roots t\ 1 . . . , t n of ip n in (—1, 1), by using 
the algorithms of Sections [5.21 15.31 respectively. We use 10 • n Gaussian nodes to evaluate 

W-It^T (639) 

and 



B ntj =icX n -^ n (l,t j ), (640) 

where ^(l,^) is defined via (|437p in Theorem [591 We observe that A n j and B n ,j appear on the 
right-hand side of ([43611 in Theorem [591 

Next, for each integer m = 0,1, . . . ,n — 1, we use the same Gaussian quadrature to evaluate 

1 4>n(t) ■ 4> m (t) dt 
1 *j 

In addition, for each integer m = 0, 1, . . . , n — 1, we compute 

|A m | 2 — |An| 2 

by using the algorithms of Sections 15.11 15.21 All the calculations are carried out in double precision. 



c 


n 


3 


|A„| 


A n ,j 


B n ,j 


10 


20 


13 


0.11487E-09 


-.25341E+01 


-.69171E-11 


500 


340 


226 


0.27418E-09 


-.19569E+01 


-.17690E-09 



Table 10: Illustration of Theorem\59l See Experiment 9. 



We display the results of the experiment in Figure [7] and Tables [TT1 [T^l In Figure [7[ we plot the 
function ip n (t)/(t — tj), corresponding to c = 10, n — 20, and j = 13, over the interval (—1, 1). We 
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¥ n (t) / (t-t) 

10 ■ ■ 




-1 -0.5 05 1 



Figure 7: The graph of ip n (t) / (t — tj) with c — 10, n — 20 and j = 13. Corresponds to Table \Tl\ See 
Experiment 9. 
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m 


J-l t-t, 


(t)dt 


|A m | 2 ^ m fe) 
|A m |MA„|- 







- 


18363E+01 


0.72463E- 


F00 


-.15543E-14 


1 


- 


28929E+01 


0.11416E- 


1-01 


-.53291E-14 


2 


- 


18299E+01 


0.72208E- 


(-00 


-.51070E-14 


3 





73457E+00 


-.28987E^ 


-00 


0.12212E-14 


4 





19270EH 


-01 


-.76041E-I 


-00 


0.37748E-14 


5 





40316EH 


-00 


-.15909E-) 


-00 


0.21094E-14 


6 


- 


14464E+01 


0.57078E- 


FOO 


-.22204E-14 


7 




10263E^ 


-01 


0.40498E- 


FOO 


0.10658E-13 


8 





11062EH 


-01 


-.43654E^ 


-00 


0.95479E-14 


9 





17030EH 


-01 


-.67204E-I 


-00 


0.99920E-14 


10 




23035E4 


-00 


0.90899E-01 


-.26645E-14 


11 




19061E4 


-01 


0.75217E+00 


-.44409E-15 


12 




91510E4 


-00 


0.36111E+00 


0.16653E-14 


13 





13774EH 


-01 


-.54355E+00 


0.11990E-13 


14 





18002EH 


-01 


-.71037E-) 


-00 


0.37748E-14 


15 




23786E-+ 


-00 


0.93863E-01 


-.97422E-14 


16 




19723E+01 


0.77830E+00 


-.11546E-13 


17 




10566E4 


-01 


0.41697E- 


FOO 


-.15987E-13 


18 





12849EH 


-01 


-.50705E-I 


-00 


-.19984E-14 


19 





19509EH 


-01 


-.76986E-I 


-00 


-.42188E-14 



Table 11: Illustration of Theorem \59\ with c = 10, n = 20 and j — 13. See Experiment 9. 



m 


rl ip n {t)ip m (t)dt 
J-l t-tj 


|A m | 2 Vm(t.,) 

|A m r-|A„|- 







0.60926E-12 


0.31814E-12 


0.33872E-15 


20 


0.43712E-01 


0.22813E-01 


0.52666E-14 


40 


-.32804E+01 


-.17120E+01 


-.84377E-13 


60 


0.85749E+00 


0.44751E+00 


-.60729E-13 


80 


-.14190E+01 


-.74055E+00 


-.91926E-13 


100 


0.84651E+00 


0.44178E+00 


-.28089E-13 


120 


0.35414E+00 


0.18482E+00 


0.47351E-13 


140 


0.53788E+00 


0.28071E+00 


-.21316E-13 


160 


-.17111E+01 


-.89302E+00 


-.35749E-13 


180 


-.93523E+00 


-.48808E+00 


0.31863E-13 


200 


-.30219E+00 


-.15771E+00 


0.48406E-13 


220 


-.51322E+00 


-.26784E+00 


0.45852E-13 


240 


-.12216E+01 


-.63753E+00 


-.11546E-13 


260 


-.10503E+01 


-.54811E+00 


-.82379E-13 


280 


0.93142E+00 


0.48609E+00 


0.84377E-14 


300 


-.55310E-02 


-.28865E-02 


0.50818E-13 


320 


0.11601E+00 


0.60544E-01 


-.13105E-12 


339 


0.14218E+01 


0.74200E+00 


-.94369E-13 



Table 12: Illustration of Theorem \59\ with c = 500. n = 340 and j = 226. See Experiment 9. 



104 



observe that this function has n — 1 roots in (—1, 1): all the roots of ip n except for tj. Obviously, 
the value of this function at tj is ip' n (tj). 

In Tables ITUl ITTl [T2l we display the results of the experiment, corresponding to c = 10, n = 20, 
j = 13 and c = 500, n = 340 and j = 226, respectively. Table [10] contains the values of the 
parameters c,n,j, as well as the quantities A n j, B n j, defined, respectively, via (|639p . (I640p above. 
Tables [TT1 1121 have the following structure. The first column contains the parameter m (an integer 
between and n— 1). The second column contains (|641j) (the left-hand side of (|436l) in Theorem 151)1) : 
in other words, this is the inner product of ip n (t)/(t — tj) with ip m . The third column contains (|642j) 
(appears on the right-hand side of ((436])). The last column contains the absolute error E m of the 



calculation of (|642[) . defined via 

nl tp n (t) ■ j) m (t) dt _ \\ m \ 2 ■ lfj m (tj) 

t — tj |A m | 2 — | A„ | 2 



E m — 



(A n .j + Bn,j) 



(643) 



(obviously, E m would be equal to zero in exact arithmetics, due to Theorem I59[) . 

We make the following observations from Tables [TU1 [11] [T^] As expected, A n j is significantly 
larger than B n j (by a factor of order |A„| _1 ). In other words, 

rl 1p n {t) ■ 1p m (t) dt _ \X m \ 2 ■ 1p m (tj) f 1 1pn(t) dt 



t - ti 



|A m | 2 -|A„|^ 



t-U 



(l + 0(|A n |)) 



(644) 



(see Theorem [55] and (jrJ59")) . l|MD|) . (EH]), (|Bl2l above). Also, in each of Tables HU [H all the 
quantities in the second and third column are roughly of the same order of magnitude (except for 
the first row in Table [T2]) . We also observe that the numerical evaluations of the left-hand side and 
the right-hand side of (|436[) in Theorem IS"9l agree up to an absolute error of order w 10 -14 . 



c 


n 


|A„| 


|-Pn,n-2| 


1000 


670 


0.93659E-11 


0.49177E-03 


1000 


690 


0.73056E-18 


0.43907E-03 


1000 


710 


0.15947E-25 


0.40076E-03 



Table 13: Illustration of Theorem \61l Corresponds to Figured 



Experiment 10. In this experiment, we illustrate Theorem IfJTI in Section 14.4.21 We proceed as 
follows. We choose, more or less arbitrarily, the band limit c and the prolate index n. Then, 
we evaluate A„, using the algorithm of Section I5~2l (in double precision). Next, for each integer 
< m < n — l,we evaluate P n , m , defined via (|448[) in Theorem l60l in Section 14.4.21 by using 
the algorithms of Sections 15.11 15.31 fin double precision). We observe that, due to Corollary [4] in 
Section 14.4.21 it suffices to consider only even values of m (since P n , m = if m is odd) . 

We display the results of the experiment in Tables [T31 HH and in Figure [S] In Figure [S] we plot 
l-fn.ml as a function of even integer m on the logarithmic scale for c = 1000 and three choices of 
n, namely, n = 670 (pluses), n = 690 (circles), and n = 710 (triangles). The value 2c/ n is marked 
with a red dashed line. In Table we display the quantities |A„| and |P n! „_ 2 |, corresponding to 
Figure [5] 

We make the following observations from Figure [S] Table [T3] and some additional numerical 
experiments. First, |P n , m | < |A„| for all m < 2c/ n (obviously, in Figure [8] we see this phenomenon 
only for n = 670, since the calculations are carried out in double precision; for n = 690, 710 and 
m < 2c/ir, \P n ,m\ < 10~ 15 ). On the other hand, for even 2c/ir < m < n, we observe that |P n ,m| 
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n, m 




800 



Figure 8: Plot of P„, m pE8")) with c = 1000 and n = 670 (crosses), n = 690 (circles), n = 710 
(triangles). The value m — Icjix is marked with a dashed line. See Experiment 10. 
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c 


n 


|An| 


max P„ >m | 

0<m<n 


c • max P n , m 

0<m<n 


50 


47 


0.26917E-07 


0.81444E-02 


0.40722E+00 


50 


53 


0.72096E-11 


0.72290E-02 


0.36145E+00 


50 


57 


0.19830E-13 


0.67353E-02 


0.33676E+00 


100 


81 


0.19431E-07 


0.48065E-02 


0.48065E+00 


100 


87 


0.18068E-10 


0.44412E-02 


0.44412E+00 


100 


93 


0.10185E-13 


0.41418E-02 


0.41418E+00 


250 


179 


0.18854E-07 


0.22730E-02 


0.56825E+00 


250 


186 


0.22556E-10 


0.14014E-02 


0.35035E+00 


250 


193 


0.17851E-13 


0.20475E-02 


0.51188E+00 


500 


339 


0.40938E-07 


0.12600E-02 


0.63000E+00 


500 


348 


0.20575E-10 


0.85073E-03 


0.42537E+00 


500 


355 


0.39965E-13 


0.11550E-02 


0.57751E+00 


1000 


659 


0.38241E-07 


0.68143E-03 


0.68143E+00 


1000 


668 


0.44256E-10 


0.49838E-03 


0.49838E+00 


1000 


677 


0.35933E-13 


0.63339E-03 


0.63339E+00 


2000 


1297 


0.41740E-07 


0.36453E-03 


0.72906E+00 


2000 


1307 


0.47570E-10 


0.35192E-03 


0.70385E+00 


2000 


1317 


0.39064E-13 


0.34212E-03 


0.68424E+00 


4000 


2572 


0.33682E-07 


0.16247E-03 


0.64987E+00 


4000 


2583 


0.37417E-10 


0.18703E-03 


0.74813E+00 


4000 


2594 


0.30728E-13 


0.14902E-03 


0.59608E+00 



Table 14: Illustration of Theorem\61l See Experiment 10. 
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grows roughly exponentially with m, reaching its maximum at m = n — 2. This maximum is 
approximately 5 • 1CP 4 , for all the three values of n (see Table [T5|) . However, Theorem IfJTl asserts 
that, for all m < n, 



I Pfl , m | — 



(645) 



c 



In other words, Theorem 1611 overestimates |-P n , m | by a factor of order n 2 . 

In Table we display some additional results of this experiment. This table has the following 
structure. The first and second column contain, respectively, the band limit c and the prolate index 
n. The third column contains |A n |. The fourth column contains 



The last column contains the value (|646l) . multiplied by the band limit c (i.e. the product of the 
first and fourth columns). 

We make the following observations from Table [TJ] and some additional experiments. First, for 
each of the seven values of c, the three indices n were chosen in such a way that |A n | is between 
10 -14 and 10 -7 . Even though c varies between 50 (the first three rows) and 4000 (the last three 
rows), the values in the last column are roughly of the same order, for all the choices of c and n. 
Moreover, these values are always between 0.3 and 0.75. This observation seems to indicate that 
Theorem loTI overestimates this quantity by 0(n 2 ) (see also (I645P and Figure 

Additional observations seem to indicate that the maximum in (|646[) is always attained at the 
largest even m between zero and n — 1 (as in Figure [8]). Also, for this value of m, all the summands 



in (|448p have been observed to have the same sign for all j = 1, . . . ,n. Thus, the inaccuracy of 
the bound in Theorem 1611 is due to overestimation of the summands (|647[) . rather than due to 
cancellation of summands with opposite signs. 

6.2 Performance of the Quadrature 

In this subsection, we report the results of numerical experiments illustrating the performance of 
the quadrature, defined in Definition [2j and whose properties are studied in Section 14.41 

6.2.1 Quadrature Error and its Relation to |A ra | 

Experiment 11. In this experiment, we illustrate Theorem IBTil in Section T4.4.2I We proceed as 
follows. We choose, more or less arbitrarily, the band limit c and the prolate index n. We evaluate X n 
as well as the nodes t±, . . . , t n and the weights Wi, ■ ■ ■ , W n of the quadrature, defined in Definition [5] 
in Section 14.41 To do so, we use the algorithms of Sections 15.21 15.31 15.41 respectively (in double 
precision). 

Then, we choose an even integer < m < n, and evaluate A m , ip m (0) and ip m (tj), for all 
j = 1, 2, . . . , n, by using the algorithms of Sections |5.2[ 15.11 (in double precision). Next, we evaluate 
P n . m , defined via (|448|) in TheoremloTil (see Experiment 10 in Section l6.1.4p . Also, we compute ||/||oo 
(see (|450l) in Theorem [501 and (|638[) in Experiment 8 in Section 16.1. 3D . Finally, we evaluate 



max{|P„ jrra | : < m < n 



!}• 



(646) 



^m(*j)^n,j(l) 



(647) 




(648) 
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m 




S m 


I V'm (1 - E <Pj) 


Cm 


error 


n 

U 


n 7f)fifi9F+00 


706R9F+00 

U . 1 UUUc/X^H^UU 


90856F-1 6 

U.Z/UOtJUiJ XU 


TQ^QQF-I 9 

U. 1 uO&iJi—i Li 


- 55^1 1 F-1 5 

. U U U X X I_J X U 


2 


49581 F+OO 


49581 F+00 


77098F-1 5 

U. 1 I Ut/OJ-J -LU 


99496F-1 

U . Zj C/t: Zj ULj X U 


- 8881 8F-1 5 

.OOOIOX/ xu 


4 


49581 F+00 


49581 F+00 

U.4:Zj OO J_ J—/ I UU 


9790nF-1 5 

U.i7 1 ZjUUXJ -LU 


37098F-10 

U.O 1 Uc"OX/ xu 


- 93870F-14 

.ZOO 1 ULv L^ 


U 


38597F+00 


38597F+00 


- 83346F-1 5 

.0iJ04:ULJ -LU 


.UlO X ULj X U 


- 1 3393F-14 

. X OUZj OLj L^ 


8 


35695E+00 


35695E+00 

U . OO U *JO±-J T UU 


- 10918E-14 


- 41 671 F-1 


- 99990E-1 5 


i n 


3351 fiF+00 

U. OUU 1UU T 


3351 6F+00 

U. OUU 1UU T UU 


U.ZjUUUOXJ -LU 


97596F-1 1 

U . cf I U ZjUXJ X X 


- 1 7908F-14 

.X I ZiUOLj LI 


12 


31 730F+00 


31 730F+00 

U.Ol / JUU T uu 


- 9^i500F-14 

. Z< U UUU-L-J Xt: 


- 97396F-1 

.C7 1 OZjUXJ xu 


111 P19F-1 

U . X X X UZjXJ X u 


14 


30901 F+00 


30901 F+00 

U.OU^UliJ T uu 


- 3^496F-14 


. X OUZj X X_J U !_7 


1 3878F-1 4 

U.XUO 1 OXJ L'i 


1 fi 


98844F+00 

U . iOO^T:i-J i^UU 


98844F+00 

U . I/OOt:t:Lv [ UU 


- 90470F-14 


- 78198F-10 

. J OXZjOXJ XU 


- 1 6653F-1 ^ 

.lUUUdX/ XU 


18 


0.27604E+00 


0.27604E+00 


-.28733E-13 


-.10967E-08 


0.42188E-14 


20 


0.26435E+00 


0.26435E+00 


-.14073E-12 


-.53714E-08 


0.90483E-14 


22 


0.25299E+00 


0.25299E+00 


0.26178E-11 


0.99913E-07 


0.94924E-14 


24 


0.24150E+00 


0.24150E+00 


0.15530E-10 


0.59274E-06 


-.66613E-15 


26 


0.22919E+00 


0.22919E+00 


-.17315E-09 


-.66085E-05 


-.72997E-14 


28 


0.21377E+00 


0.21377E+00 


-.53359E-09 


-.20365E-04 


0.14710E-14 


30 


0.18075E+00 


0.18075E+00 


0.55489E-08 


0.21178E-03 


-.51903E-14 


32 


0.10038E+00 


0.10038E+00 


-.62071E-08 


-.23690E-03 


-.70915E-14 


34 


0.27988E-01 


0.27988E-01 


-.88231E-07 


-.33675E-02 


0.10113E-13 


36 


0.49822E-02 


0.49818E-02 


0.40165E-06 


0.15330E-01 


0.29751E-14 


38 


0.70503E-03 


0.70008E-03 


0.49503E-05 


0.18894E+00 


-.13444E-13 



Table 15: Illustration of the proof of Theorem \60\ with c = 50 and n = 40. See Experiment 11. 



where the functions ipi, . . . , ip n are those of Definition [2] in Section l4~4l 

We display the results of the experiment in Table [15] The data in this table correspond to 
c = 50 and n = 40. Table [T5l has the following structure. The first column contains the even integer 
parameter m, which varies between and n — 2. The second column contains A m ')/' m (0) (we observe 
that 

A m ^ m (0)= / i> m {t)dt, (649) 



due to (I57|) in Section l2~Tj) . The third column contains the quantity S m , defined via the formula 

IA„J 2 



|A m | 2 -|A„|= 



ic\ n Pn,m + ^ ^m(tj) ■ W j } . (650) 
.7 = 1 



The fourth column contains the integral (|648p . The fifth column contains the number £ m , defined 
via the formula 



1 f 1 



\I\\ 



J ^ m {t)- fl- (651) 



(We observe that, due to (|452l) in Theorem [501 £ m equals to the value in the third column, divided 
by ||/||oo- The latter does not depend on m, and is equal to 0.26201E-04, for c = 50 and n = 40.) 
The last column contains the difference between the value in the third column and the sum of the 
values in the fourth in fifth columns (due to the combination of (|648|) . (|649[) . (|650l) and (|452[) in 
Theorem 1601 this quantity would be zero in exact arithmetics). 
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We make the following observations from Table [T5l First, 

S m = J Vn.(i)-fo>i(t) + ••• + ¥>*(*))*, (652) 

due to the combination of (|650|) and Theorem IBTfl We observe that S m is close to X m ip m (0) for small 
to, but coincides with the latter only in two digits for m = 38. 

Second, we observe that the value in the fourth column (see (16481) ) is grows from w 10~ 16 at 
to = up to « 5 • 1(T 6 at m = 38. 

Third, we observe that, due to the combination of (|651[) and Theorem IBTil 

1 f 1 

U = 7Tj7— / I(t) ■ il> n (t) ■ i> m (t) dt, (653) 
IM Hoc J -i 

where / is that of Theorem I5B1 in Section l4.3.3l Theoretically, |£ m | is bounded from above by 1 (see 
the proof of Theorem l60|) . However, in fact, |£ m | is significantly smaller than one for small values of 
to, though £ m w 0.2 for to = 38. 

Next, the value in the last column, that would be zero in exact arithmetics, serves as a test of 
the accuracy of the calculation. We observe that this value is of order 10 -14 , for all to. Finally, we 
note that A m - m (O) is always positive and monotonically decreases with to. 



TO 




/ i>m - S Wjtpm 


(tj) 




Cn.m/l-^nl 





0.70669E+00 


-.44409E-15 







26389E- 


■04 





20432E+00 


2 


0.49581E+00 


-.16653E-15 







26333E- 


04 





20389E+00 


4 


0.42581E+00 


-.13323E-14 







26314E- 


04 





20375E+00 


6 


0.38527E+00 


-.21649E-14 







26303E- 


04 





20366E+00 


8 


0.35695E+00 


-.22760E-14 







26296E- 


04 





20361E+00 


10 


0.33516E+00 


-.16653E-14 







26290E- 


04 





20356E+00 


12 


0.31730E+00 


-.23870E-14 







26285E- 


04 





20352E+00 


14 


0.30201E+00 


-.24980E-14 







26281E- 


04 





20349E+00 


16 


0.28844E+00 


0.11102E-14 







26277E- 


04 





20346E+00 


18 


0.27604E+00 


-.59230E-13 







26274E- 


04 





20344E+00 


20 


0.26435E+00 


0.83716E-12 







26271E- 


04 





20342E+00 


22 


0.25299E+00 


-.89038E-11 







26268E- 


04 





20339E+00 


24 


0.24150E+00 


0.76862E-10 







26265E- 


04 





20337E+00 


26 


0.22919E+00 


-.65870E-09 







26262E- 


04 





20335E+00 


28 


0.21377E+00 


0.45239E-08 







26253E- 


04 





20327E+00 


30 


0.18075E+00 


-.19826E-07 







26282E- 


04 





20350E+00 


32 


0.10038E+00 


0.68548E-07 







26276E- 


04 





20345E+00 


34 


0.27988E-01 


-.33810E-06 







26849E- 


04 





20789E+00 


36 


0.49822E-02 


0.27232E-05 







28516E- 


04 





22080E+00 


38 


0.70503E-03 


-.22754E-04 







72700E- 


04 





56291E+00 



Table 16: Illustration of Theorem \60\ with c = 50 and n = 40. See Experiment 12. 



Experiment 12. In this experiment, we illustrate Theorems IBTil 1621 in Section T4.4. 21 We proceed 
as follows. We choose, more or less arbitrarily, band limit c and prolate index n. We evaluate A„, 
as well as the nodes ti,...,t n and the weights Wi, . . . , W n of the quadrature, defined in Definition [5] 
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in Section POl To do so, we use, respectively, the algorithms of Sections [5.11 15.21 15.31 15.41 (in double 
precision). Then, we choose an even integer < m < n, and evaluate A m , ^ m (0), and ip m (tj) for all 
j = 1, . . . , n, using the algorithms of Sections 15.21 15.11 (in double precision). 

We display the results of this experiment in Table [TBI The data in this table correspond to c = 50 
and n = 40 (the same as for Table [TBI in Experiment 11). Table [TBI has the following structure. The 
first column contains the even integer m, that varies between and n — 2. The second column 
contains A m ^m(0). The third column contains the difference 

n 

A m V™(0) - ' Wj. (654) 

J'=l 

The fourth column contains the number C„. m , defined via the formula 

|A "'^ ||/||oo + |A„|- fJMl^m(0)|+c|P„, m |V (655) 



Am| 2 / Vl-^ml 

where ||/||oo and P n , m are defined, respectively, via (|450[) and (|448D in Theorem [SD] (see also Exper- 
iment 8 in Section [6. 1 .31 and Experiment 10 in Section [6.1 .41) . Note that (|655[) is the right-hand side 
of (|449[) in Theorem l60l The fifth column contains C n . m /\\ n \. 

We make the following observations from Table [Jj)J We note that (|654[) in the third column is 
the error of the quadrature rule of Definition [21 used to integrate tp m over (—1, 1) (see also ([3"T]) in 
Section l2~Tj) . The absolute value of this error is close to the machine precision for small m, and 
grows up to PS 2 • 10~ 5 for m = 38. For all values of m, the absolute value of (|654[) is bounded by 
C n ,m (the fourth column), in agreement with Theorem 1601 We also observe that C„ im is of the same 
order of magnitude for all values of m (as opposed to (|654l) ). Moreover, C n . m is always smaller than 
|A„| (in this case, |A„| = 0.12915E-03). More specifically, C n , m is between 0.2 • |A„| and 0.6 • |A„|, 
for all the values of m (see the last column) . 

The behavior of the quadrature error (|654[) in the third column is explained with the help 
Experiment 10 and Table [TBI in Experiment 11, as follows. Due to (|453[) in the proof of Theorem IBTil 

n 

A m <M0)-^Vmfe)-Wj = 

|2 Vm(0) , A |A„| 2 



IK]' ■ +11- ll^lloo + ic\ n Pn, m , (656) 

A m V \Am\ J 

where ^ m is defined via (|65ip in Experiment 11. The first summand in the right-hand side of (|656p 



grows as m increases. The behavior of the second summand in the right-hand side of (|656[) depends 
on £ m , which is close to zero for small values of m and close to one for large values of m (see the 
fifth column in Table [TBI . Finally, the last summand in the right-hand side of (|656[) is also expected 
to grow with m (compare to Figure [5] in Experiment 10). 

To conclude, C„ jm , defined via (|655[) . significantly overestimates the quadrature error (I654p for 



small values of m. On the other hand, when m is close to n, C njTn is a fairly tight bound on (|654[) 
In Theorem 1621 we provide an upper bound on C n ^ m (and hence on the quadrature error (|654[T 
which is independent on m, namely, 



A m V>m(0) - ^ 1pm(tj) ■ Wj 



< c n 



<\\ n \- ( 24 • log ( -L) +6-x n ). MioT ! 



|A„| 
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However, the logarithmic term in (|657p is due to the inaccuracy of Theorem [521 in Section l4.3.2l (see 
Experiment 7 in Section l6.1.3[) . Also, the term 6 • \n m (1657[) is due to the inaccuracy of Theorem loTI 
in Section [4.4.21 (see Experiment 10 in Section |6.1.4[) . In other words, numerical experiments seem 
to suggest that the quadrature error (|654[) is bounded by |A„|, for all even < to < n. 



Quadrature error 
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Figure 9: The quadrature error \J ip m — X^rofe) ■ Wj\ as a function of even to < n, for four 
different values of n and c = 10000. See Experiment 12. 



In Figure [9l we display the results of the same experiment with different choice of parameters c 
and n. Namely, we choose c = 10000 and plot \ m ip m {fy as a function of even < m < 6425, on the 
logarithmic scale (solid blue line). In addition, we plot the absolute value of the quadrature error 
(|654p . as a function of to, for four different values of n: n — 6393 (red dashed line), n = 6401 (red 
circles), n = 6414 (red triangles), and n — 6425 (red pluses). The corresponding values of |A n | are 
displayed in Table [TTJ 



n 


6393 


6401 


6414 


6425 


|An| 


0.43299E-07 


0.54119E-09 


0.33602E-12 


0.52616E-15 



Table 17: Values of \X n \ for c — 10000 and different choices of n. 



We make the following observations from Figure [9l First, A m i/' m (0) is approximately a constant 
for to < 2c/7r, and decays roughly exponentially with to for to > Icjix. Also, for each value of n, 
the quadrature error (I654[) is essentially zero for to < 2c/ir, and its absolute value increases roughly 
exponentially with to for to > 2c/ V. Nevertheless, the absolute error of the quadrature error is 
always bounded from above by |A n |, for each n. See also Tables [TBI IT51 and Conjectured] below. 

We strengthen the observations above by repeating this experiment with several other values of 
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c 


n 


m 


AmV'm(O) 




|A„| 


250 


179 


178 


0.28699E-07 


-.52496E-08 


0.18854E-07 


250 


184 


182 


0.68573E-09 


-.38341E-10 


0.16130E-09 


250 


188 


186 


0.14108E-10 


-.68758E-12 


0.30500E-11 


500 


339 


338 


0.52368E-07 


-.13473E-07 


0.40938E-07 


500 


345 


344 


0.37412E-09 


-.86136E-10 


0.27418E-09 


500 


350 


348 


0.12148E-10 


-.99816E-12 


0.35537E-11 


1000 


659 


658 


0.42709E-07 


-.14354E-07 


0.38241E-07 


1000 


665 


664 


0.51665E-09 


-.15924E-09 


0.43991E-09 


1000 


671 


670 


0.52494E-11 


-.15024E-11 


0.42815E-11 


2000 


1297 


1296 


0.41418E-07 


-.17547E-07 


0.41740E-07 


2000 


1304 


1302 


0.77185E-09 


-.15036E-09 


0.37721E-09 


2000 


1311 


1310 


0.31078E-11 


-.11386E-11 


0.28754E-11 


4000 


2572 


2570 


0.54840E-07 


-.15493E-07 


0.33682E-07 


4000 


2579 


2578 


0.43032E-09 


-.20771E-09 


0.46141E-09 


4000 


2587 


2586 


0.28193E-11 


-.12805E-11 


0.29164E-11 


8000 


5119 


5118 


0.43268E-07 


-.26751E-07 


0.52899E-07 


8000 


5128 


5126 


0.50230E-09 


-.16395E-09 


0.33442E-09 


8000 


5136 


5134 


0.50508E-11 


-.15448E-11 


0.32132E-11 


16000 


10213 


10212 


0.42725E-07 


-.30880E-07 


0.56568E-07 


16000 


10222 


10220 


0.69663E-09 


-.28201E-09 


0.52821E-09 


16000 


10231 


10230 


0.34472E-11 


-.22162E-11 


0.42902E-11 



Table 18: Relation between the quadrature error and |A„|. See Experiment 12. 
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band limit c and prolate index n. The results are displayed in Table [T5] This table has the following 
structure. The first and second column contain, respectively, the band limit c and the prolate index 
n. The third column contains the even integer < m < n (the values of m were chose to be close 
to n). The fourth column contains X m tp m (0). The fifth column contains the quadrature error (I654p . 
The last column contains |A„|. 

We make the following observations from Table [T5] First, for each of the seven values of c, the 
three indices n were chosen in such a way that |A n | is between 10~ 12 and 10 -7 . The values of the 
band limit c vary between 250 (the first three rows) and 16000 (the last three rows). For each n, 
the value of m is chosen to be the largest even integer below n. This choice of m yields the largest 
A m VVre(0) and the largest quadrature error (|654l) among all m < n (see also Table [TB)). Obviously, 
|A m | and |A„| are of the same order of magnitude, for this choice of m. We also observe that, for all 
the values of c, n, m, the absolute error of the quadrature error (|654j) is bounded from above by |A n | 
(and is roughly equal to |A„|/2). In other words, the upper bound on the quadrature error, provided 
by Theorem 1621 (see (|657|1 ). is somewhat overcautious. 

We summarize these observations in the following conjecture. 

Conjecture 2. Suppose that c > is a positive real number, and n > 2c/7r is an integer. Sup- 
pose also that < m < n is an integer. Suppose furthermore that t\, . . . ,t n and Wi, . . . , W n are, 
respectively, the nodes and weights of the quadrature, introduced in Definition^ in Section \4-4\ Then, 



Remark 26. Conjecture^ provides a stronger inequality than that of Theorem \6 6 A On the other 
hand, Conjecture has been only supported by numerical evidence, while Theorem \62\ has been 
rigorously proven. 

Experiment 13. In this experiment, we demonstrate the performance of the quadrature, intro- 
duced in Definition [5] in Section POI on exponential functions. We proceed as follows. We choose, 
more or less arbitrarily, the band limit c and the prolate index n. We evaluate the quadrature 
nodes t\,...,t n and the quadrature weights W\, . . . , W n , by using, respectively, the algorithms of 
Sections [5.31 15.41 (in double precision). Also, we evaluate |A„|, by using the algorithm in Section I5T21 
(in double precision). Then, we choose a real number a > 0, and evaluate the integral of e %cax over 
— 1 < x < 1 via the formula 




(658) 




(659) 



Also, we compute an approximation to (|659l) . by evaluating the sum 



n 




(660) 



i=i 



Finally, we evaluate the error of this approximation, that is, 



2 sin(ac) 



1 1 




(661) 
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e icxa _ £ e icxt. w 




0.2 0.4 0.6 0.8 1 

Figure 10: The quadrature error (|66ip with c = 1000, n = 650. See Experiment 13. 



e icxa _ £ e icxt. w 




0.5 1 1.5 2 

a 

Figure 11: The quadrature error (|661[) with c = 1000, n = 650. See Experiment 13. 
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In Figures [TU1 QTJ wc display the results of this experiment. The band limit and the prolate 
index were chosen to be, respectively, c = 1000 and n — 650. This choice yields A„ = -.21224E-04. 
In these figure, we plot the quadrature error (|661[) as a function of the real parameter a. Figure [TU] 
corresponds to < a < 1, while Figure ITTI corresponds to < a < 2. 

We make the following observations from Figures [10J[TTJ For < a < 1, the absolute value of the 
quadrature error (|661|) is bounded by 4 • 10~ 9 ~ 10 • |A„| 2 . The largest quadrature error is obtained 
when a is close to 1. On the other hand, for 1 < a < 2, the absolute value of the quadrature error 
(|66ip is significantly larger, and is of order |A„|. The largest quadrature error is obtained when a is 
close to 1. 

These observations admit the following (somewhat imprecise) explanation. Suppose that a > 
is a real number. Due to ([37|l and Theorem [1] in Section |2~T1 



e <«» = J2 *mi>m(a)Tp m (x), (662) 
m=0 

for all real — 1 < x < 1 (we note that while e lacx is not a bandlimited function of —1 < x < 1, it 
does belong to L 2 [—1, 1]). Moreover, 

f e — dx = = g A 2 > m (a)^ m (0). (663) 



We combine (|56T|) . p)3"l) . to obtain 

2 sin(ac) 



• cos(icat,) = 



oo / n \ 

^ A m ^m(a) A m V,n(0) - ^ W^^mfe) ■ (664) 

m=0 y j=l / 

We recall (see Experiment 12) that, for small values of m, the quadrature error (|654[) is very small 
compared to \X n \- On the other hand, for those values of to < n that are close to n, the quadrature 
error (|654[) is of order |A n |. Therefore, roughly speaking, 

71—1 I n \ 

\ m iP m {a) A m Vm(0) -J]l^i(tj) = O (|A n | 2 • Vn-i(a)) . (665) 

m=0 \ j=l / 

On the other hand, due to the fast decay of |A m |, we expect 

oo / n \ 

£ \ m ipm{a) X m ip m {0) - W^ m {t ) = O (|A„| 2 • ^„(o)) . (666) 

m=n y j = l J 

If < a < 1, then \^ n (a)\ = 0{^/n) (see Theorems IT21 [HI [T4l in Section |2~T1) . We combine this 
observation with (|665l) . f|666[) to conclude that the quadrature error (|661l) is expected to be of the 
order |A„| 2 • y/n. 

If, on the other hand, 1 < a < 2, then |i/> n (a)| = O (|A„| _1 ) (see, for example, Theorem IMl in Sec- 
tion SXU Theorem g3] in Section rC2~2l Theorem gS] in Section 143^21 Experiment 1 in Section |6~1~T1 
Experiment 6 in Section r6.1.2j) . We combine this observation with (|665[) . (|666p to conclude that, in 
this case, the quadrature error (|661l) is expected to be of the order |A„|. 

We summarize this crude analysis, supported by the observations above, in the following conjec- 
ture. 
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Conjecture 3. Suppose that c > is a real number, and that n > 2c/tt is an integer. Suppose also 
that t\ , . . . , t n and W± W n are, respectively, the nodes and weights of the quadrature, introduced 
in Definition^ in Section \4-4\ Suppose furthermore that — 1 < a < 1 is a real number. Then, 



p ic«i 



^e~~> ■W J =0(\\ n \ 2 -Vn~). (667) 



Experiment 14. In this experiment, we illustrate Theorems IB"4l 1651 in Section T4.4. 31 We proceed 
as follows. We choose, more or less arbitrarily, the band limit c > and the accuracy parameter 
e > 0. Then, we use the algorithm of Section 15721 to find the minimal integer m such that |A m | < e. 
In other words, we define the integer n\(e) via the formula 

ni(s) = mm{m > : |A m |<e}. (668) 

Also, we find the minimal integer such that the corresponding bound on the quadrature error, 
established in Theorem R)2l in Section 14.4.21 is less that e (see also (|657[) in Experiment 12). In other 
words, we defined 77.2(e) via the formula 

n 2 (e) = min jm > : |A m | • (^24 • log (j^j) + 6 • Xm^j < £ j • (669) 

Then, we define the integer 713(e) via the formula (|493[) in Theorem l64l In other words, 

/2c a(e) , / 16ec \\ ,„„„s 
n 3( e) = floor [- + -£■ log (^y, j J (670) 

where a(e) is defined via (|492[) in Theorem[64] Finally, we define the integer 714(e) via the right-hand 
side of (|506p in Theorem 1551 In other words, 

n 4 (e) = floor ^ + ^10 + | - log(c) + ^ • log ^ • log (|) J . (671) 

In both (|670p and (|67ip . floor(a) denotes the integer part of a real number a. 

We display the results of this experiment in Table Q1J1 This table has the following structure. 
The first column contains the band limit c. The second column contains the accuracy parameter e. 
The third column contains 711(e), defined via (|668D . The fourth column contains 712(e), defined via 
(|669[> . The fifth column contains n 3 (e), defined via ([6701) . The sixth column contains 714(e), defined 
via (1671[) . The seventh column contains |A ni ( e )|. The last column contains |A n2 ( e \|. 

Suppose that c > is a band limit, and n > is an integer. We define the real number Q(c,n) 
via the formula 



Q(c, n) — max 



ip m (t) dt-^2 tPrn(tj) ■ Wj 
3=1 



0<m<n-l} , (672) 



where t\, . . . , t n and W\, . . . , W n are, respectively, the nodes and the weights of the quadrature, 
defined in Definition [2] in Section 14.41 In other words, this quadrature rule integrates the first n 
PSWFs up to an error at most Q(c, n). 

We make the following observations from Table IT9l We observe that Q(c, 77.1(e)) < e, due to the 
combination of Conjecture [2] in Section f6. 2.11 and (|668[) . (I672p . In other words, numerical evidence 
suggests that the quadrature of order 77.1(e) will integrate the first ni(e) PSWFs up to an error at 
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25 
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OOZ 
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U. 04:0Z c7 JJj Z(J 
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OUU 


io- 


50 


4oo 


AAK 
440 


OU I 


7 A O 
1 4z 


U.ozoy I-Cj-DU 


fl QQOI 7T? 
U.oozl ( Hj-00 


1UUU 


io- 


10 


add. 
000 


bo < 


oUo 


oo4 


U.yooozrj-lU 


n noo/i 7T7 1 7 
U.yzy4 (ej-LI 


1000 


io- 


25 


707 


725 


875 


942 


0.97844E-25 


0.14241E-31 


1000 


io- 


50 


767 


783 


981 


1120 


0.39772E-50 


0.56698E-57 


2000 


io- 


10 


1305 


1330 


1467 


1500 


0.95177E-10 


0.25349E-17 


2000 


io- 


25 


1351 


1373 


1550 


1619 


0.86694E-25 


0.27321E-32 


2000 


io- 


50 


1418 


1438 


1675 


1818 


0.88841E-50 


0.22795E-57 


4000 


io- 


10 


2581 


2610 


2768 


2804 


0.70386E-10 


0.64396E-18 


4000 


io- 


25 


2632 


2658 


2862 


2935 


0.57213E-25 


0.53827E-33 


4000 


io- 


50 


2707 


2730 


3007 


3154 


0.56712E-50 


0.88819E-58 


8000 


io- 


10 


5130 


5163 


5344 


5383 


0.59447E-10 


0.22821E-18 


8000 


io- 


25 


5185 


5216 


5450 


5526 


0.87242E-25 


0.16237E-33 


8000 


io- 


50 


5268 


5296 


5614 


5765 


0.95784E-50 


0.23927E-58 


16000 


io- 


10 


10225 


10264 


10468 


10509 


0.63183E-10 


0.37516E-19 


16000 


io- 


25 


10285 


10321 


10585 


10664 


0.85910E-25 


0.41416E-34 


16000 


io- 


50 


10377 


10409 


10769 


10923 


0.51912E-50 


0.56250E-59 


32000 


io- 


10 


20413 


20457 


20686 


20730 


0.62113E-10 


0.12818E-19 


32000 


io- 


25 


20478 


20519 


20815 


20897 


0.78699E-25 


0.12197E-34 


32000 


io- 


50 


20577 


20615 


21018 


21176 


0.96802E-50 


0.15816E-59 


64000 


io- 


10 


40786 


40837 


41092 


41139 


0.89344E-10 


0.28169E-20 


64000 


io- 


25 


40857 


40903 


41232 


41318 


0.66605E-25 


0.39212E-35 


64000 


io- 


50 


40964 


41008 


41454 


41616 


0.85451E-50 


0.28036E-60 



Table 19: Illustration of Theorems \ &4\ [Eh\ See Experiment 14. 
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most e (see Remark |2T)]) . On the other hand, we combine Theorem 1521 in Section [4.4.21 with (I669[) . 
(16721) . to conclude that the quadrature of order /12(e) has been rigorously proven to integrate the 
first 712(e) PSWFs up to an error at most e. In both Theorem [52] and Conjecture we establish 
upper bounds on Q(c,n) in terms of |A n |. The ratio of |A ni ( e )| to |A„ 2 ( e )| is quite large: from about 
10 6 for c = 250 and e = lCT 10 , 1(T 25 , 1(T 50 (see the first three rows in Table US), to about 10 10 for 
c = 64000 and and e = 10~ 10 , 10~ 25 , 10" 50 (see the last three rows in Table [15]). On the other hand, 
the difference between 712(e) and ni(e) is fairly small; for example, for e = 10~ 50 , this difference 
varies from 10 for c = 250 to 23 for c = 4000, to merely 44 for as large c as c = 64000. 

As opposed to 771(e) and 712(e), the integer 77.3(e), defined via (|670[) . is computed via an explicit 
formula that depends only on c and e (rather than on |A n | and Xn> that need to be evaluated 
numerically). This formula is derived in Theorem 1641 bv combining Theorem 1621 with some explicit 
bounds on |A„| and \n in terms of c and 77. The convenience of (|670[) vs. (|668[) . (|669p comes at a 
price: for example, for e = 10 -50 , the difference between 723(e) and 712(e) is 123 for c = 250, and 446 
for c = 64000. However, the difference 773(e) — 772(e) is rather small compared to c: for example, for 
e = 10~ 50 , this difference is roughly 4 • (log(c)) , for all the values of c in Table IT9"1 

Furthermore, we observe that 714(e) is also computed via an explicit formula that depends only 
on c and e (see ([671]) ) . This formula is a simplification of that for 773(e), derived in Theorem [65] 
Thus, not surprisingly, 774(e) is greater than 773(e), for all the values of c and e. 

We summarize these observations as follows. Suppose that the band limit c and the accuracy 
parameter e > are given. In Theorem l641 we prove that 77 > 773 (e) implies that the quadrature error 
Q(c, 77), defined via (|672[) . will be at most e (for the quadrature of order 77, defined in Definition [5] 
in Section 23]) . On the other hand, numerical evidence suggests that Q(n,c) < e also for all the 
values of n between 771(e) and 773(e) (see Experiment 12). In this experiment, we observed that the 
difference between 773(e) and ni(e) is relatively small compared to c (roughly of order (log(c)) ). 



6.2.2 Quadrature Weights 

Experiment 15. In this experiment, we illustrate the results of Section 14.4.41 fin particular, The- 
orem [S7] Corollary [5] and Remark [14]) . We proceed as follows. We choose, more or less arbitrarily, 
band limit c and prolate index 77. Then, we compute the quadrature nodes t\, . . .,t n as well as 
tfj' n (ti), . . . , ip' n (t n ) , by using the algorithm of Section [531 We evaluate tp' n (0), using the algorithm 
of Section 15.11 Next, we evaluate the quadrature weights Wi, . . . , W n , by using the algorithm of 
Section [54] Also, for each j = 1, . . . , 77, we evaluate the sum 

4 ( 673 ) 



where Qk{t) is the fcth Legendre function of the second kind, defined in Section l2~2l and a£ is the 
fcth coefficient of the Legendre expansion of ip n , defined via ((54"]) in Section |2~2"1 fsee Theorem 1571 and 
Section [54]) . To evaluate (|673[) numerically, we use only 2N first summands, where N is an integer 
of order 77 (see (|587[) in Section I5"4"j) . All the calculations are carried out in double precision. 

We display the results of this experiment Table [201 The data in this table correspond to c = 40 
and 77 = 41. Table[20]has the following structure. The first column contains the weight index j, that 
varies between 1 and 21 = (77 + l)/2. The second column contains Wj. The third column contains 
the difference between Wj and (|673[) . The last column contains the difference 

W j ^iSfpl m 

(see Remark [T4"|) . 
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3 




W j + 2-$ n (t j )/TP' n (t j ) 


Wj- 


(^(i,)) 2 -(l-t?) 


1 





7602931556894E-02 


0.00000E+00 




55796E-11 


2 





1716167229714E-01 


0.00000E+00 




55504E-10 


3 





2563684665002E-01 


0.00000E+00 




21825E-12 


4 





3278512460580E-01 


0.00000E+00 




11959E-09 


5 





3863462966166E-01 


0.16653E-15 





82238E-11 


6 





4334940472363E-01 


0.22204E-15 




16247E-09 


7 





4713107235981E-01 


0.22204E-15 





11270E-10 


8 





5016785516291E-01 


0.19429E-15 




18720E-09 


9 





5261660773966E-01 


0.26368E-15 





10495E-10 


10 





5460119701692E-01 


0.29837E-15 




20097E-09 


11 





5621699326080E-01 


0.17347E-15 





81464E-11 


12 





5753664411864E-01 


0.12490E-15 




20866E-09 


13 





5861531690539E-01 


0.10408E-15 





55098E-11 


14 





5949490764741E-01 


0.23592E-15 




21301E-09 


15 





6020725336886E-01 


0.13184E-15 





31869E-11 


16 





6077650804037E-01 


0.18041E-15 




21545E-09 


17 





6122088420703E-01 


0.48572E-16 





14361E-11 


18 





6155390478472E-01 


0.83267E-16 




21675E-09 


19 





6178529976346E-01 


0.11102E-15 





36146E-12 


20 





6192162112196E-01 


0.48572E-16 




21732E-09 


21 





6196665001384E-01 


0.00000E+00 


0.00000E+00 



Table 20: Quadrature weights (|33"4]) with c = 40, n = 41. A„ = i0.69857E-08. See Experiment 15. 
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In Figure [T2l we plot the weights Wj, displayed in the second column ol Table |2"01 For j > 21, 
the weights are computed via symmetry considerations. Each Wj is plotted as a red dot above the 
corresponding node tj. 

We make the following observations from Table [30] First, all the weights are positive (see 
Theorem [751 and Remark [111]) . Moreover, Wj grow monotonically as j increases to (n + l)/2. Also, 
due to the combination of Theorems (57J EH1 in Section 14.4.41 the value in the third column would be 
zero in exact arithmetics. We observe that, indeed, this value is zero up to the machine precision, 
which confirms the correctness of the algorithm of Section [5^41 (We note that, for j = 1, 2, 3, 4 and 
j = 21, this algorithm, in fact, does evaluate Wj via (|673[> . and hence this value in the corresponding 
rows is exactly zero). Finally, we observe that, for all j, the value (|674[) in the last column is of the 
order |A„|, in correspondence with Remark IT4"1 

Quadrature weights 
0.08, ■ ■ ■ i 



0.06 



0.04 



0.02 







-1 -0.5 




0.5 



Figure 12: The quadrature weights W\, . . . , W n with c = 40, n = 41. See Experiment 15. 
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